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Abstract 



The increasing number of discovered extra-solar planets opens a new opportunity for 
studies of the formation of planetary systems. Their diversity keeps challenging the 
long-standing theories which were based on data primarily from our own solar system. 
Resonant planetary systems are of particular interest because their dynamical config- 
uration provides constraints on the otherwise unobservable formation and migration 
phase. 

In this thesis, formation scenarios for the planetary systems HD128311 and HD45364 
are presented. N-body simulations of two planets and two dimensional hydrodynami- 
cal simulations of proto-planetary discs are used to realistically model the convergent 
migration phase and the capture into resonance. The results indicate that the proto- 
planetary disc initially has a larger surface density than previously thought. 

Proto-planets are exposed to stochastic forces, generated by density fluctuations in a 
turbulent disc. A generic model of both a single planet, and two planets in mean motion 
resonance, being stochastically forced is presented and applied to the system GJ876. 
It turns out that GJ876 is stable for reasonable strengths of the stochastic forces, but 
systems with lighter planets can get disrupted. Even if a resonance is not completely 
disrupted, stochastic forces create characteristic, observable libration patterns. 

As a further application, the stochastic migration of small bodies in Saturn's rings 
is studied. Analytic predictions of collisional and gravitational interactions of a moon- 
let with ring particles are compared to realistic three dimensional collisional N-body 
simulations with up to a million particles. Estimates of both the migration rate and 
the eccentricity evolution of embedded moonlets are confirmed. The random walk of 
the moonlet is fast enough to be directly observable by the Cassini spacecraft. 

Turbulence in the proto-stellar disc also plays an important role during the early 
phases of the planet formation process. In the core accretion model, small, metre- 
sized particles are getting concentrated in pressure maxima and will eventually un- 
dergo a rapid gravitational collapse to form a gravitationally bound planetesimal. Due 
to the large separation of scales, this process is very hard to model numerically. A 
scaled method is presented, that allows for the correct treatment of self-gravity for a 
marginally collisional system by taking into account the relevant small scale processes. 
Interestingly, this system is dynamically very similar to Saturn's rings. 
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Chapter 



Introduction 



And if the fixed Stars are the Centers of other hke 
systems, these, being form'd by the hke wise counsel, 
must be all subject to the dominion of One, [...]. 

Isaac Newton, General Scholium, translated by Motte, 1729 

In 1713 Isaac Newton wrote this sentence in an essay attached to the third edition of his 
famous Principia Mathematica. In other words, he expected to see planets around other 
stars. Almost three centuries later, astronomers discovered the first planet beyond our 
own solar system, a so-called exo-planet. The number of known planets has increased 
rapid ly ever since. To date, 461 extra-solar planets have bee n discovered ( Schneidej 



20101) ■ At least 10% of all nearby solar type stars host planets ( ICumming et al.i,2008.) . 
With this tremendous observational success, it is now the theoreticians' turn to explain 
the discovered systems. One important aspect is to find out if, and if so why, these 
systems formed differently compared to the solar system. 

In this chapter, first the discoveries of exo-planets in recent years are presented. 
Then, suggested formation scenarios of planets, planetary systems and their evolution 
are reviewed. 

This thesis discusses stochastic phenomena in a range of astrophysical systems. 
The analytic model presented in chapter |3] is the key in understanding the effects of 
stochastic forces and turbulence in those systems. It forms the basis of the physical 
understanding in many celestial systems in which stochastic forces are present. 

A detailed formation scenario of the planetary system HD45364 is presented in 
chapter |2j In chapter [3] another formation scenario is presented, this time for the 
planetary system HD128311. Both systems are resonant systems and their dynamical 
states provide important constraints on their formation history. 

IID45364 formed most likely in a massive disc and had a phase of rapid convergent 
migration. On the other hand, the current observed orbital parameters of the system 
HD128311 are consistent with the formation in a strongly turbulent disc. 

In chapter |U these formalism developed in chapter [3] is applied to Saturn's rings 
and a moonlet. Saturn's rings also exert stochastic forces. The rings, together with 
embedded moons, resemble a small scale version of the proto-planetary disc. 

Finally, the issue of numerical convergence in simulations of planetesimal formation 
is discussed in section [51 Planetesimals are likely to form in turbulent proto-stellar discs 
via gravitational instability. It turns out that the system can be simulated consistently 
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only if the relevant small scale processes are included. The dynamical evolution is then 
very similar to Saturn's rings, except that the final clump is gravitationally bound. 

In chapter [6] we summarise the results. The main numerical codes that have been 
used in chapters 121 El IH and El are described in the appendices [E] and 

1.1 Methods of detecting extra-solar planets 

For thousands of years, the observation of planets was limited to the major planets of 
our own solar system. Since then, our knowledge of the solar system has been growing 
constantly and has reached an overwhelming magnitude. However, the fundamental 
question of whether we inhabit a special place in the universe or not hasn't been an- 
swered yet. This is deeply linked to our spiritual desire to know if there are lifeforms 
on other planets. Luckily, we are living in a time of great technological and scientific 
progress and it is not unreasonable to assume that those questions can be answered 
within the next 50 years. 

The first steps have already been taken, namely the discovery of planets beyond 
our own solar system. Various groups around the world have successfully detected exo- 
planets using different techniques. Each method has both advantages and disadvantages 
which will be summarised in this section. This is important because the sparse infor- 
mation that we get from these observations determines the predictability of theoretical 
studies. 

Radial velocity measurements. Most planets have been discovered by the radial 
velocity (RV) method. A periodic Doppler shift in the spectral lines of the host 
star can be measured with high precision spectroscopy. The Doppler effect occurs 
because the star is moving around the centre of mass (of both the star and the 
planet). Only the radial part of this movement is measurable. The period of 
the oscillation is simply the planet's orbital period. The mass of the star can be 
calculated by stellar evolution models. If we assume that we observe the system in 
the plane of the planet's orbit, we can then use Newton's third law MqVq = nip Vp, 
where Mq, rup, Vq and Vp are the masses and velocities of the host star and the 
planet respectively, to calculate the planet's mass. 

This method biases the detection of massive planets on short orbits (hot Jupiters) 
as the gravitational influence of the planet on the star (and therefore the measured 
Vq) is bigger in those cases. Note that there is a degeneracy in the inclination i 
because it is not possible to measure the non-radial part of the velocity. Thus, we 
can only get a lower limit on the mass of the planet and the mass might actually 
be larger: 



Because of this degeneracy, only a limited number of parameters can be obtained 
by the RV method. Finding orbital parameters of multiplanetary systems be- 
yond the planet masses and the orbital periods is even more challenging. Each 
new planet adds 7 degrees of freedom, such that solutions are highly degener- 
ate. Especially the fltting of the eccentricities is very unrehable (see discussion in 
chapter [2]). Furthermore, it is difficult to flnd precise orbital parameters for very 
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long period planets as it can take many years to sample a full period of the light 
curve. 

T he first extra-solar plan et discovered by the radial velocity method is 51 Pegasi 
b flMavor fc QTldo3ll995f ). 



Transit light curves. Another method that has been very successful is the transit 
method. This is mainly due to space based missions such as CoRoT and Kepler. 
The idea is to find planets by observing variations in the star light caused by 
transits of planets in front of the host stars. This effect also occurs in our own 
solar system, known as Mercury and Venus transits. However, transits of extra- 
solar planets are rare because the host star, the exo-planet and the Earth must be 
aligned exactly in one line for a transit to occur. The first exo -planet discovered 
by the transit method is OGLE-TR-56 b (IKonacki et al.ll2003l ). 

The transit method is capable of measuring the density p of exo-planets, a pa- 
rameter inaccessible to the RV method. It is even possible to do spectroscopy on 
the atmosphere of the planet during the transit and the secondary transit (when 
the planet is occulted by the star) to create a temperature map of the exo-planet 
flKnutson et al.llioOTh . 

Furthermore, the Rossiter-McLaughlin effect allows one to measure the sky pro- 
jected angle between the orbit and the rotation axis of the host star. To do that, 
one has to combine transit and RV measurements. A recently submitted paper 
suggests that a large number of planets might be on retrograde orbits (Triaud et 
al, in preparation). 

Dedicated ground and space based missions promise to detect a large number of 
planets in the near future. Furthermore, by observing tiny variations in the transit 
timing (TTV) and transit duration (TDV), it might be possible to find other 
planets in systems in which only one pl anet is transiting. E ven the possibility to 
discover exo-moons has been discussed (IKipping et al.l 120091 ). 



Gravitational microlensing. Planets can also be discovered by gravitational mi- 
crolensing. If the star is aligned in one line with the Earth and a bright object in 
the background (e.g. another star), the light from the background object bends 
around the star and can be detected by an increase in luminosity. A planet in 
an orbit around the star disturbs the light curve and theoretical models can de- 
termine the planet's mass and some orbital parameters. This event happens only 
once per star, so good timing and a global collaboration is needed to perform a 
continuous measurement of the light curve. The orbital parameters cannot be 
determined with high precision because one only obtains a snapshot of the sys- 
tem without any dynamical evolutio n. The first exo-planet was discovered by 



gravitational microlensing in 2004 by iBond et al. 



Direct imaging. In 2004 IChauvin et al.l reported the first detection of a giant planet 
candidate by direct imaging. Since then, several planetary systems have been 
imaged. These are all massive planets, far away from the host star (hundreds of 
AU). Maybe the most interesti ng of those planet s is Fomalhaut b, which seems to 
be embedded in a debris disc (IKalas et al.ll2008l ). As all those planets have long 
orbital periods, it is, similar to the gravitational microlensing method, difficult to 
determine the precise orbital configuration, especially in multi-planetary systems. 
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Pulsar timing. Despite the recent success of the detection methods p resented above, 



the fi rst exo-planet was discovered using the pulsar timing method (jWolszczan fc Frail 



19921 ) by measuring slight variations in the regular timings from a pulsar. This is 
sensitive down to very small mass planets. However, pulsars are rare compared 
to normal stars and the focus has shifted away from the pulsar timing method. 

1.2 Observed objects 

Extrasolar planets have been observed around a variety of parent stars from pulsars 



to solar- type stars to M-dwarfs (see e.g. IChauvin et al.l 12004 IWolszczan fc Frail! Il992l ) 



indicating that planet formation is common and successful in a broad range of envi- 
ronments. Almost all extra-solar planetary systems are distinct from the solar system. 
Many of the detected objects are so-called hot Jupiters. Their size and mass is com- 
parable to Jupiter, but their orbits are very close to their host star 0.05 AU). Most 
methods described above bias the detection in favour of these objects. 

The existence of hot Jupiters was very surprising because in the solar system all 
gas giants are located beyond several AU. This still results in difficulties for planet 
formation theories, although planetary migration is one solution, as described below. 

Other observed objects are very heavy and are more likely to be brown dwarfs rather 
than planets. The International Astronomical Union (lAU) defines an extra-solar planet 
as follows: 

Objects with true masses below the limiting mass for thermonuclear fusion of 
deuterium (currently calculated to be 13 Jupiter masses for objects of solar 
metallicity) that orbit stars or stellar remnants are planets (no matter how 
they formed). The minimum mass/size required for an extra-solar object to 
be considered a planet should be the same as that used in our Solar System. 
flWGESPll2003[ ) 



Thus, many of the massive objects, at separations of hundreds of AU, detected by direct 
imaging, could actually be brown dwarfs. 

In cases where the planet's mean density can be observed, evolutionary models seem 
to be in broad agreement with observations, although many hot Jupiters have a large 
surface temperature and are subject to tidal heating. Their atmospheres are not well 
understood and first studies show a broad range of possible configuratio ns. Observations 



furthermore suggest that their densities might be rather low (see e.g. [Anderson et al 



20101) 



The Holy Grail for exo-planet hunters is another Earth-like planet that orbits the 
host star within the habitable zone. At the present day, the planet that is closest to 
one Earth mass is Gliese 581 e which has a mass of approximately 2 Earth masses 



flMavor et all 120091 ). 



1.3 Formation of planets 
1.3.1 Accretion disc 

Stars form out of giant gas clouds that become gravitationally unstable. Because of 
angular momentum conservation, an accretion disc forms around every new star. As 
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the name suggests, an accretion disc will eventually accrete most of its material onto 
the star, leaving only a small fraction of it as planets or a debris disc in orbit around 
the star. Observations of accretion rates in pro to-planetary discs suggest that accretion 



timescales are of the order of a few 10^ years (lEnoch et al.l 12008 



Usually one assumes that the Navier-Stokes equations are a good approximation to 
describe the fluid motion within the disc. However, to reproduce the measured accretion 
rates, a purely molecular viscosity is not effici ent enough. I t is generally believed that 
this large viscosity originates from turbulence (jPringldll98ll ). The mass fiux in a steady 
state Keplerian disc is then given by 



1.2) 



where z/ is the effective kinematic viscosity and S the surface density. From observations 
of proto-stars, a value of z/ ~ 10~^AU^/ (yr/27r) has been inferred, but with considerable 
uncertainty. 



The magneto rotational instability (MRI, iBalbus &: Hawlevlll99lf ) is the most likely 
candidate to be responsible for the an omalous value of u which is often characterised 
using the IShakura Sz Svunvaev a parameter, defined as a = z//(c^/f2), with 

being the local sound speed and 2tt/VL being the orbital period. For the most likely 
situation of small or zero net magnetic fiux, recent MHD simulations have indicated 
10^2 - 10^1 This value is very uncertain when realistic values of the actual 



a 



tra nsport coefficients are employed due in no small pa rt to numerical resolution issues 



see 



Fromang fc Papaloizoull2007l : iFromang et al.ll2007l ) . Which parts of proto-planetary 



discs are adequa tely ionised, or constitute a dead zone, is also an issue (jGammidll996 ; 
Sano et al.ll200oh . 

Putting aside all those issues for a moment, the MRI always creates density fiuctu- 
ations in the disc, resulting in a stochastic force that is exerted on embedded planets 
which are otherwise decoupled from the gas motion (ignoring migration, which is a lam- 
inar effect and acts on timescales much longer than one orbit). In this thesis, we do not 
attempt to simulate the MRI directly and rather describe it in an empirical way. We 
therefore avoid all problems mentioned above and can understand the physical scaling 
of our results. Only two quantities are needed for our analytic model, the root mean 
square val ue of the stochastic g ravitational force and the corresponding auto correla- 
tion time. iLaughlin et al.l ( 120041 ) propose a more complicated model, in which random 
modes with random decay times create a gravitational potential that is supposed to 
mimic the MRI. Our, much simpler and more intuitive description, allows us to survey 
a large parameter space and understand the resulting physical processes, as discussed 
in chapter [31 



1.3.2 Minimum mass solar nebula and snowline 



A standard model of a proto-planetary disc assumes a steady state and a vertical 
equilibrium (see also appendix IE.1.2|) . The gas and dust components are furthermore 
assumed to be fully mixed. The resulting disc is fiared, although still geometrically 
thin. For sufficiently low mass accretion rates (< 10~ ^Mf7i/yr) the dominant heating 
so urce is ste l lar ir radiation ( IChiang fc GoldreichI 119971 ). 



Havashil ( 1l98lh prescribes the so-called minimum mass solar nebula (MMSN) which 



is comprised of just enough mass to make all planets of the solar system. This model 
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became the standard disc model and has been used excessively in recent years with 
different normalisations. One can describe the surface density and temperature profiles 
as 



Sn I - 

^0 



T(r) 



ni- 



1.3) 



Typical normalisation values at tq = lAU are Sq = 1700 gem an d Tf 



280 K. 



Standard values of p range from to 5/3, those of q from 1/2 to 3/4 ( lHavashilll981 



Cuzzi et al.lll993l ). The typical mass of the disc is about one percent of the stellar mass. 
Although the choice of disc model is essential for m ost aspects of planet form ation, little 
emphasis has been placed on alternatives (see e.g. iDeschI 120071 : ICridall2009l ). 

To form planets a nd plane tary cores, dust is an important ingredient (see also be- 
low). In the model of Hayashil (q = 1/2), the radius at which temperatures drop below 
170 K is around 2.7 AU. At larger radii, the temperatures are low enough for water ice 
to exist. This transition radius is called the snow line. Recent studies, including more 
detailed models of accretiona l heating and radiative transfer show that the snow line 
could come as close as 1 AU fISasselov fc Lecai]l2000f ). 



1.3.3 Planet formation mechanism 

Planets are believed to form in proto-stellar discs as a natural by-product of star forma- 
tion. These discs are made of the same material as the star itself: gas and dust. Current 
theories give two possible explanations of the formation of giant planets inside the disc. 
Both models favour planet formation at large radii, beyond the snowline. However, it 
is important to keep in mind that the process of planet formation itself is not directly 
observable (yet), leaving theory and numerical simulations to fill in the blanks between 
observations of hot circumstellar discs around young stars and planets orbiting main 
sequence stars. 



Core accretion 



One of the most important unanswered questions in the theory of planet formation 
is to find out what the mechanism for planetesimal formation is, i.e., the process by 
which the building blocks of planets are formed. In the core accretion model, solid 
components of the disc stick toget her, forming bigger and bigger objects until a core o f 
about 15 Earth masses is formed (IBodenheimer fc Pollacklll986l : iPollack et al.lll996al ). 
The proto-planet can then start to accrete gas and form a gas giant. 

Again, there are two main t heories for planetesi mal formation via the core accretion 
model: mut ual collisions fe.g.. iHavashi et al.l 119771 ) and gravitational instability in the 
dust layer (iGoldreich fc WardI Il973l ). In the first hypothesis, dust particles grow as 
the result of accretion-dominated collisions. Although the formation of planetesimals 
by mutual collisions is consistent with meteoritic evidence, the collision speed between 
dust particles (or aggregates) must be much slower than the t ypical velocity disper sion 
in a standard proto-stellar disc to avoid destructive collisions (IBlum fc Wurnj|2008l ). In 
addition, the planetesimal formation process is so slow that metre-sized particles are in 
danger of spiralling into the star before growing large enough to decouple from the gas 
( IWeidenschillingI 1 19 77bl ) . Even if km-sized planetesimals were able to form, they would 
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be in danger of being ground down again by mutual collisions (IStewart &: Leinhardt 

2nn9h . 

Gravitational instability is often considered to be a solution to most of these prob- 
lems because the intermediate sizes are avoided all together. In this theory, the dust 
layer becomes dense enough for the Keplerian shear and velocity dispersion of the 
dust particles to be unable to support the dust against its own self gravity. The dust 
then collapses into clumps that eventually cool via drag forces and mutual collisions 
into planetesimals. Many different groups have been working on this subject. Until 
recently, the focus has been on quiet, non-turbulent, and low density regions of the 



accretion disc (see e.g., iMichikoshi et al.ll2009l . 120071 iTanga et al.l 120041 ). 

However, the turbulent gas in the proto-planetary nebula stirs the dust, which 
increases the velocity dispersion of the dust particles. Several ideas have been proposed 
to overcome the turb ulence-induced mi xing of the dust particles and create localised 
clumps. For example, ICuzzi et al.l (120081 ) suggest that the same turbulence that stirs the 
dust on larger s cales may also collect the dust particles on small scales. A similar idea 
was proposed by lJohansen et al.l (120071 ) . in which dust partic les are localised into clump s 
promoted by both turbulence and the streaming instability (lYoudin fc Goodmarul2005l ). 
These dense clumps then become gravitationally unstable. A third hj^othesis suggests 
that large structures, such as vo rtices, may be able to collect and protect du st particles 
from the turbulent background (IBarge fc Sommerialll995l : iLyra et al.ll2009l ). 

In chapter \5\ we focus on the gravitational collapse in a very dense and turbulent 
region of the proto-planetary disc. We look carefully at the numerical requirements 
of modelling gravitational instability accurately and test the validity of using super- 
particles in a high density region. The results suggest that some of the early results of 
graviational collapse from other authors may have been too optimistic. 



Gravitational fragmentation 

In this m odel, gas p lanets form directly as a result of gravitational instability within 
the disc (jBosd l200ll ): no solid core is needed. However, the planet can accrete dust 
particles later on, and hence form a core. 

For the gravitationa l instability to occur, it is necessary to have a Toomre Q pa- 
rameter of order unity (jToomrdll964l ). The precise criterion depends strongly on the 
thermodynamics on the disc, especially the cooling time. However, it is unlikely that 
the physical conditions in a proto-planetary disc are compatible with this constraint. 
If so, this is most likely to occur at large radii (~ 100 AU). It therefore cannot be the 
formation mechanism for most discovered exo-planets. 



1.4 Evolution of planetary systems 



1.4.1 Migration 

Many observed planets are very close to their host star. It is implausible that they 
have formed at such small radii (see section [1.3.21) . even taking into account the strong 
selection effect of discovering close in exo-planets (see section [TTTl) . 

Tidal interactions with the proto-planetary disc give the planet some radial mobi l- 
ity and might therefore be the solution to this problem (IGoldreich fc Tremaindll980l ). 
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The mutual angular momentum exchange depends on many parameters of both the 
planet and the disc. Several population synthesis simulatio ns were able to r eproduce 
the observed mass-period distribution using migration (e.g. Ilda fc LinI l2004l ). Planet 
migration can be classified in the following four main categories. 



Type I 

Planets are completely embedded in a proto-planetary disc for sufficiently low mass. 
Density waves are excited in the disc, both interior and exterior to the planet, at so- 
called Lindblad resonances. The waves carry angul ar momentum and produce a. torque 
on the planet that leads to planetary migration ( iGoldreich fc Tremaind Il979l ). The 
direction of migration depends on parameters of the disc model such as the gradients 
of the surface density, sound speed and scale height (I WardI 119971 ). In most cases, the 
outer torque is bigger than the inner one and the planet migrates inward. This effect 
is called type I migration. 

If the migration rate is very fast, low mass planets are in danger of spiralling into 
the star within the disc lifetime. The precise speed and direction in a realistic disc 
model are still subject to debate, where recently the focus h as been on non-isothermal 
equations of state and the effect s of co-rotation torques (jPaardekooper et al.l 12010 : 
Paardekooper fc Papaloizoull2009l ). 



Type II 

If the planet is massive enough, it can open a gap in the disc. For most disc models 
about a Jupiter mass is needed to clear a clean gap. The process is similar to the gaps 
opening in Saturn's rings with moons orbiting within the gap. The gap establishes 
a flow bar rier to the disc rnateria l, effectively locking the planet to the viscous disc 
evolution ( iLin fc Papaloizoul 119861 ). Thus, the planet still migrates, the regime being 
called type II migration. However, the migration timescale is set by the disc evolution 
timescale which is in general much longer than the type I migration timescale. 



Type III 

In a disc with a relatively large surface density (several times the MMSN), the density 
distribution in the co-orbital region of the planet can be asymmetric. This leads to a 
self- sustained torque that is proportional to the migration rate. In this regime, called 
type III migration, the planet can fall inward on a timescale much shorter than the 
disc evolution time obtained for type II migration, even shorter than the timescales 
associated with type I migration. The planet moves faster than the disc can respond 
to its perturbation, thus maii itaining an asymmetry . The net torque can cause both 
inward or outward migration f Peplihski et al. 2008a , b, 3). 

In chapter [2] we present the flrst observable indication that type III migration is 
indeed responsible for shaping planetary systems. 



Type IV 

The different migration regimes have all been studied in quiet non-turbulent discs. 
However, a proto-planetary disc is thought to be at least partially turbulent (see section 
ll.3.ip . Assuming the simplest scenario in which the effects of turbulence and the net 
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migration are separable, one can describe the stochastic migration as a diffusion process 
in a distribution of planets on top of the net orbital migration, described above. We 
call this regime type IV. 

This idea leads to a Fokker-Planck description (lAdams fc BlochI 120091 ). There are 
two limits to this approach. The first, in which the net migration dominates and the 
effects of turbulence are negligible (mostly for heavy planets). The second, in which the 
net migration can be ignored and the turbulence driven diffusion dominates on short 
timescales (mostly for small mass planets and planetesimals) . We are in an intermediate 
regime if the migration and diffusion timescales are comparable. It is unlikely that the 
interplay of stochastic migration and net migration plays an important role in determin- 
ing the final semi major axis of planets because the associated timescales would have to 
be very similar. However, resonant systems are more sensitive to stochastic forces and 
might therefore provide observational constraints on the strength of turbulence that 
was present at early times. 

So far, only preliminary work has been done on studying th e effect of the turbulence 
on resonances (lAdams et al.l l2008l : iNelson fc Papaloizoul |200J) . If the turbulence is ac- 
tive long enough, it will eventually kick the planets out of resonance. This might happen 
even if the stochastic forces are not strong enough to cause significant orbital migration. 
However, the timescale needed for destroying a resonance is not well constrained. New 
results on this issue are presented in chapter [31 



1.5 Mean motion resonances 



Resonances in celestial m echanics can occur when there is a simple numerical relation 
between two frequencies ( Murray &: Dermott 2000l ). For example, a planet could feel 
a periodic gravitational force from another planet that is a multiple of its own orbital 
frequency. This can either lead to an unstable situation in which angular momentum is 
exchanged until the resonance ceases to exist or a planet gets ejected from the system, 
or a stable configuration. 

Resonances can form whe n dissipatiye forc es act on the planets, for example in a 
proto-planetary disc (see e.g. iMalhotral (119931 ) and appendix [B]) . When there are two 
(or more) planets in the disc, the migration rates might differ for various reasons. One 
possibility is that the planets are in different migration regimes (see sect ion [1.4. II) . The 
resulting differential migration changes the ratio of orbital periods Pi and the ratio of 
semi major axes aj. These ratios are related by Kepler's third law, (Ti/T2)^ = (01/02)^. 

The planets can become locked into a mean motion resonance (MMR) if the ratio 
of orbital periods gets close to a rational number q/p with small integers q and p, e.g. 
1/2, 1/3 or 2/3. The numbers p and q can be interpreted as the number of completed 
orbits after the same time. Some of these resonances are stable and planets can stay in 
resonance over a long period of time, migrating together and keeping the ratio of their 
orbital periods constant. Beside the presence of Hot Jupiters, resonance capturing is 
admitted as strong evidence that planets undergo a phase of orbital migration. This 
idea was su ccessful in explaini i ig the observed resonan t multi-planet systems GJ876 and 
55 Cancri ( Lee fc Pealel[200l[ : [Snell grove et al.[[200l[ l which are in a 2:1 resonance In 
chapter [21 we discuss the first successful formation scenario for a system that is in a 3:2 
resonance. 
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Chapter 



The dynamical origin of the 
multi-planetary system HD45364 

Truth is ever to be found in simplicity, and not in the 
multiphcity and confusion of things. 

Issac Newton, unpublished manuscript, Prank E. Manuel 1974 



The recently discovered planetary system HD45364, which consists of Jupiter- and 
Saturn-mass planets, is very likely in a 3:2 mean motion resonance. The standard 
scenario for forming planetary commensurabilities involves convergent migration of two 
planets embedded in a proto-planetary disc. However, when the planets are initially 
separated by a period ratio larger than two, convergent migration will most likely lead 
to a stable 2:1 resonance, incompatible with current observations. 

Rapid type III migration of the outer planet crossing the 2:1 resonance is one possible 
way around this problem. Here, we investigate this idea in detail. We present an 
estimate of the required convergent migration rate in section 12.1.21 and confirm this 
with N-body simulations in section 12.21 and hydrodynamical simulations in section 12.31 
If the dynamical history of the planetary system had a phase of rapid inward migration 
that forms a resonant configuration as we suggest here, then we predict that the orbital 
parameters of the two planets will always be very similar and should show evidence of 
that. 

We use the orbital parameters from our simulation to calculate a radial velocity 
curve and compare it to observations in section 12.51 Our model provides a fit that is 
as good as the previously reported one. However, the eccentricities of both planets are 
considerably smaller and the libration pattern is different. Within a few years, it will be 
possible to observe the planet-planet interaction directly and thus distinguish between 
these different dynamical states. 

This is the first prediction of orbital parameters for a specific extra-solar planetary 
system derived from planet migration theory alone. It provides strong evidence on how 
the system formed. 
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2.1 The planetary system HD45364 



2.1.1 The standard formation scenario and its problems 



The planets in this system have masses of mi = 0.1906Mjup and m2 = 0.6891Mjup and 
are orbiting the star at a distance of ai = 0.6813 AU and 02 = 0.8972 AU, respectively 
( ICorreia et al.l 120081 ). The period ratio is close to 1.5. This alone does not imply that 
the planets are in a 3:2 mean motion resonance. However, a stability analysis shows 
that the three body system (two planets and a star) is only stable if the planets are in a 
3:2 mean motion resonance. Furthermore, the region of greatest stability also contains 



the b est statistical, Keplerian fit to the radial velocity measurements (ICorreia et al 
2008h . 



In the core accretion model, a solid core is fir stly formed by dust aggregation. After 
a critical core mass is at tained (i MizunJ '1980^), the proto-planet accretes a gaseous 
envelope from the nebula f lBodenheim er & Pollack 1986). 

As explained in more detail in section ll.3.2[ the planets have most likely formed 
further out in cooler regions of the proto-stellar disc as water ice, which is an important 
ingredient for dust aggregation, can only exist bey ond the snow line. The snow line is 
generally assumed to be at radii larger than 2 AU ( ISasselov &: Lecarll2000l ). 

When the planets have obtained a substantial part of their mass, they migrate due 
to planet disc interactions (see section 11.4. ip . Although the details of this process are 
still being hotly debated, the existence of many resonant multiplanetary systems and 
hot Jupiters supports this idea and suggests that planets preferably migrate inwards. 
Both planets in the HD45364 system are inside the snow line, implying that they should 
have migrated inwards significantly. 

The migration rate depends on many parameters of the disc such as surface density, 
viscosity, and the mass of the planets. The planets are therefore in general expected 
to have different migration rates, which leads to the possibility of convergent migra- 
tion. In this process the planets approach o rbital commensu rabilities. If they do this 
slowly enough, resonance capture may occur (jGoldreiclj|l965l ) . after which they migrate 
together maintaining a constant period ratio thereafter. 

Studies by several authors have shown that when two planets, either of equal mass or 
with the outer one more massive, undergo differential convergent migration, capture into 
a mean motion cor nmensurability is expe cted provided that the convergent migration 
rate is not too fast (ISnellgrove et al.ll200ll ) . The observed inner and outer planet masses 
are such that, if (as is commonly assumed for multiplanetary systems of this kind) the 
planets are initially separated widely enough that their period ratio exceeds 2, a 2:1 
commensurability is ex pected to form at low migration rates (e.g. iNelson &: Papaloizou 
2OO2I; lKleY_eraLfl2004 | ). 



Pierens fc Nelson! fl2008h have studied a similar scenario where the goal was to re- 
semble the 3:2 resonance between Jupiter and Saturn in the early solar system. They 
also find that the 2:1 resonance forms in early stages; however, in their case the inner 
planet had the higher mass, whereas the planetary system that we are considering has 
the heavier planet outside. In their situation the 2:1 resonance is unstable, enabling the 
formation of a 3:2 resonanc e later on, and the migration rate can stall or even reverse 
flMasset fc Snelkrovel EoOlh . 

For the planetary system HD45364, this standard picture poses a new problem. 
Assuming that the planets have formed far apart and were not much smaller during 
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the migration phase, the outcome is almost always a 2:1 mean motion resonance, not 
3:2 as observed. The 2:1 resonance that forms is found to be extremely stable. One 
possible way around this is a very rapid convergent migration phase that passes quickly 
through the 2:1 resonance. 



2.1.2 Avoiding the 2:1 mean motion resonance 

We found that, if two planets with masses of the observed system are in a 2:1 mean 
motion resonance, which has been formed via convergent migration, this resonance is 
very stable. An important constraint arises, because at the slowest migration rates, 
the 2:1 resonance is expected to form rather than the observed 3:2 commensurability 
provided the planets start migrating outside any low-order commensurability. 

We can estimate the critical relative migration timescale Ta^crit above which a 2:1 
commensurability forms with the condition that the planets spend at least one libration 
period migrating through the resonance. The resonance's semi-major axis width Aa 
associated with the 2:1 resonance can be estimated from the condition that two thirds 
of the mean motion difference across Aa be equal in magnitude to 27r over the libration 
period. This gives 



Aa 



uJifa2 
n2 



(2.1) 



where a2 and n2 are the semi major axis and the mean motion of the outer planet, 
respectively. The l ibration period 2ti I uju can be expresse d in terms of the orbital 
parameters (see e.g. lGoldreichlll965l : iRein fc Papaloizoull2009l . and also chapter [3]) but is 
here measured numerically, for convenience. If we assume the semi-major axes of the two 
planets evolve on constant (but different) timescales |ai/di| = r^,! and |a2/a2| = ''"a,2, 
the condition that the resonance width is not crossed within a libration period gives 



l/7-a,l - 1/Ta,; 



> 27r 



a2 



Wi/Aa 



27r 



n2 



(2.2) 



If 



to pass through the 2:1 MMR. 

If the planets of the HD45364 system are placed in a 2:1 resonance with the inner 
planet located at 1 AU, the libration period 27r/a;// is found to be approximately 75 yrs. 
Thus, a relative migration timescale shorter than r^^crit ~ 810 yrs is needed to pass 
through the 2:1 resonance. For example, if we assume that the inner planet migrates 
on a timescale of 2000 years, then the outer planet has to migrate with a timescale of 



ra,2,crit < 576 yrs. 



(2.3) 



2.2 N-body simulations 

We ran A^-body simulations to explore the large parameter space and confirm the 
estimate from the previous section. In an N-body simulation all objects (stars, planets) 
are treated as point masses interacting only gravitationally with each other. To model 
the planet-disc interaction, one can explicitly add dissipative and stochastic forces. 
Thus, the total force acting on an object i is a sum of the following terms 

Fi Fi gravity ~l~ Fi indirect ~l~ Fi a-damping ~l~ Fi c-damping ~l~ Fi stochastic l^'^) 
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The first term is due to the gravitational interaction. In units where G = 1 we have 



-^i gravity / ^ ^ j ~ ^TTS' (^■^) 



where we sum over all particles in the system, except the i-th particle. In the heliocentric 
coordinate system which is used here, the star is located at a fixed position at r = 0. 
That brings about an additional indirect term as we are in an accelerated, non-inertial 
frame, 

Fi indirect = " mirrij-^. (2.6) 

The next two terms in eqation 12.41 are due to the laminar planet-disc interaction. In 
a first approximation the interaction damps the semi-major axis a and the eccentricity 
e on time scales and Tg, respectively. The exact terms depend on the orbital pa- 
rameters of the planet and are given in appendix [Bl It is common practice to define 
the ratio between the timescales K = Ta/Tf.. To obtain those timescales, we have to 
compare the N-body simulations to full hydrodynamical simulations which will be done 
in the following section. An N-body simulation is much faster than a hydrodynamical 
simulation and thus can be integrated over a longer time span. This is essential be- 
cause resonance capturing and migration act on timescales ~ 1000 — 10000 years and 
we might even want to integrate over millions of years. The last term in equation 12.41 
adds stochastic forces, simulating the turbulent nature of proto-planetary discs. The 
implementation is discussed in chapter |3l 

Newton's second law together with equation 12.41 forms a system of ordinary differ- 
ential equations. To solve them, a new N-body code has been developed which is highly 
modular and easily expandable. It incorporates different modules for stochastic forces, 
migration, data output and time-stepping. 

The choice of integrator depends on the problem. The following algorithms have 
been implemented: 

Runge-Kutta Method (RK4) The classical, fourth order Runge-Kutta method is a 
widely-used standard integrator. It is an explicit method that needs four function 
evaluations per time-step. The main disadvantage of the RK4 method is the fixed 
time-step. We have to set a specific value at the beginning of the computation 
that is not refined later on. 

Runge-Kutta-Fehlberg Method (RKF45) The Runge-Kutta-Fehlberg r nethod is 



a fift h order explicit method with an embedded fourth order method (iFehlbere 



19691 ). We can use the two different results to estimate the numerical error and a 
new time-step. We repeat the step with a smaller time-step if the error is larger 
than a specified limit eps. 

Midpoint Method The midpoint method belongs to the class of second order Runge- 
Kutta methods. Due to the low order it is not efficient to use it on its own, but 
it is used as a sub-timestep integrator by the Bulirsch-Stoer Method. 

Bulirsch-Stoer Method (BS) This method fIStoer fc Buhrschll2002l : IPress et al.lll992[ ) 



is based on a stepwise extrapolation. For each time-step we calculate the new po- 
sitions and velocities several times with different sub-time-steps Ati using the 
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modified midpoint method. We then perform an extrapolation to a perfect sub- 
step in the limit where At — )■ 0. This allows the use of very large time-steps while 
still obtaining a high accuracy. We can also use the extrapolation to estimate the 
error and thus make the method adaptive. 



The code has been used in lRein fc Papaloizod (120091 ) and lRein et al.l (l2010l ). usually 
with the Bulirsch-Stoer integrator and a precision of e = 10~^^. A symplectic integrator 
has not been implemented, although it would have better long term conservation prop- 
erties. This is because the integrator would formally not be symplectic any more once 
velocity dependent forces have been added. This is the case in all simulations presented 
here, which include either migration or stochastic forces. Note that in some cases it is 
possible to find a symplectic integrator for sys tems with velocity dependent forces, for 
example for Hill's equations flQuinn et al.ll2010f ). Such an integrator (a modified version 
of the leap-frog algorithm) has been used in shearing sheet calculations of planetary 
rings, presented in chapter HI 



2.2.1 Parameter space survey 

With the new N-body code, we can now easily scan the parameter space and confirm the 
analytic estimate of the critical migration rate needed to capture into the 3:2 resonance. 

We place both planets on circular orbits at Oi = 1 AU and 02 = 2 AU initially. 
The migration timescale for the inner planet is fixed at 1 = 2000 yrs, while the 
migration timescale for the outer planet is varied. In figures 12. 12.21 and 12.31 we plot 
the period ratio P2/P1 and the eccentricities of both planets as a function of time t 
for different migration timescales Ta,2- In all cases, there is a sharp transition of the 
final resonant configuration from 2:1 to 3:2 at around Ta^2 ~ 550 yrs. This value agrees 
closely with the analytic estimate given by equation 12.31 Note that once the planets 
are in resonance, the eccentricities quickly reach an equilibrium value as the planets 
migrate adiabatically. 

In figure I2.2l the eccentricity damping corresponds to = 30 and is therefore three 
times stronger than in figure 12.11 A value of = 30 is rather high and probably 
unphysical. It can be seen that a smaller eccentricity damping timescale pushes the 
boundary towards a longer migration time-scale. Although this effect is rather weak 
(a three times stronger eccentricity damping shifts the boundary by only 4%), it can 
be easily understood within the toy model presented above. The eccentricities have to 
rise quickly, while the planets are within the resonance width. This is more difficult for 
short eccentricity damping timescales. 

In figure 12.31 the eccentricity damping corresponds to A' = 1 and is therefore ten 
times weaker than in figure 12.11 Again, a value of i^' = 1 is probably unphysical. This 
results in a rapid rise of eccentricities as soon as the planets are in resonance. Systems 
that are in a 3:2 resonance become unstable within a few thousand years. On the other 
hand, systems that are in a 2:1 resonance remain stable, although eccentricities are high 
(ei ~ 0.5). 

The results show clearly that, if the planets begin with a period ratio exceeding 
two, to get them into the observed 3:2 resonance, the relative migration time has to be 
shorter than what is obtained from the s tandard theory of t ype II migration applied to 



these planets in a standard model disc ( iNelson et al.ll2000l ). In that case one expects 
this timescale to be larger than 10^ years, being effectively the disc evolution timescale. 
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(a) Period ratio PijPx as a function of time (y-axis) and migration 
timescale of the outer planet Ta;i (a;-axis). 
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(b) Eccentricity of the inner planet ei as a function of time (y-axis) and 
migration timescale of the outer planet ra,2 (x-axis). 
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(c) Eccentricity of the outer planet 62 as a function of time (y-axis) and 
migration timescale of the outer planet Ta^i (a;-axis). 

Figure 2.1: Parameter space survey. The migration timescale of the inner planet is 
Ta^i = 2000 yrs. The eccentricity damping is given through K = TalTf. = 10. 
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(a) Period ratio P2IP1 as a function of time (y-axis) and migration 
timescale of the outer planet Ta^2 (a;-axis). 
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(b) Eccentricity of the inner planet ei as a function of time (y-axis) and 
migration timescale of the outer planet Ta^2 (a;-axis). 
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(c) Eccentricity of the outer planet 62 as a function of time (y-axis) and 
migration timescale of the outer planet Ta^2 (a;-axis). 

Figure 2.2: Parameter space survey. Same parameters as in figure 12.11 except K 

Ta/Te = 30. 
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(a) Period ratio P2/P1 as a function of time (y-axis) and migration 
timescale of the outer planet Ta,2 (a;- axis). 
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(b) Eccentricity of the inner planet ei as a function of time (y-axis) and 
migration timescale of the outer planet Ta^2 (a;-axis). 
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(c) Eccentricity of the outer planet 62 as a function of time (y-axis) and 
migration timescale of the outer planet 2 (x-axis). 



Figure 2.3 

Ta/Te = 1. 

the left). 



Parameter space survey. Same parameters as in figure 12.11 except K = 
Many systems eject one or more planets, indicated by the white colour (on 
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Table 2.1: Parameters for some of the hydrodynamic simulations. The second column 
gives the disc aspect ratio, the third column the softening length in units of the local 
scale height and the fourth column gives the initial surface density. The fifth and 
sixth columns give the number of radial and azimuthal grid points used and finally the 
simulation outcome is indicated in the last column. 



However, it is possible to obtain the required shorter migration ti mescales in a massive 
disc i n which the planets n iigrate in a type III regime (see e.g. iMasset &: Papaloizou 



20031 : iPeplihski et al.ll2008al ). In this regime, the surface density distribution in the co- 



orbital region is asymmetric, leading to a large torque which is able to cause the planet 
to fall inwards on a much shorter timescale than the disc evolution time obtained for 
type II migration. Other parameters such as the eccentricity damping timescale do not 
have a strong effect on capture probabilities. 

We explore the feasibility of this scenario in more detail in the following sections. 



2.3 Hydro dynamical simulations 



Two-dimensional, grid-based hydrodynamic simulations of an accretion disc with two 
gravitationally interacting planets were performed to test the rapid migration hypoth- 
esis. The simul a tions performed here are similar in concept to those performed by 
Snellgrove et al.l (l200l[ ) of the resonant coupling in the GJ876 system that may have 
been induced by orbital migration resulting from inter action with th e proto-planetary 
disc. We performed studies using the FARGO code (IMassetl l2000l ) with a modified 
locally isothermal equation of state. Those runs are indicated by the letter F. 

Willy Kley also ran simulations including viscous heating and radiative transport 
using the RH2D code . Those runs are indicated by the letter R and more details can be 
found in IPein et all fl2010h . 

Further simulations have been performed with a new hydrodynamics code called 
Prometheus, which is described in appendix [El The code is similar to FARGO, using 
operator splitting, a staggered grid and fast orbital advection. However, it it can be 
run in two dimensions (cylindrical grid) and in three dimensions (spherical coordinates). 
Furthermore, it includes a particle and an MHD module. 



2.3.1 Initial configuration and computational set up 

We use a system of units in which the unit of mass is the central mass M^,, the unit of 
distance is the initial semi-major axis of the inner planet, ri, and the unit of time is 
{GM^/r^)^^/"^ . Thus the orbital period of the initial orbit of the inner planet is 27r in 
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these dimensionless units. The parameters for some of the simulations we conducted, 
as well as their outcomes, are given in table 12.11 

In all simulations presented here, the mass ratio of the inner and outer planet is 
qi = 2.18 ■ 10~^ and q2 = 7.89 ■ 10~^, respectively. The mass ratios adopted are those 
estimated for HD45364. The initial separation of the planets is r2/ri = 1.7. The 
viscosity is chosen to be z/ = 10~^. The simulations have the inner and outer boundary 
of the 2D grid at r, = 0.25 and Tq = 3.0, respectively. In all simulations, a smoothing 
length b = 0.6 H has been adopted, where H is the disc thickness at the planet's 
position. This corresponds to about 4 zone widths in a simulation with a resolution 
of 768 X 768. The role of the softening parameter b that is used in two-dimensional 
calculations is to account for the smoothing that would result from the vertical structure 
of the disc in three dimensional calculations (e.g. Masset et al 2006). We find that the 
migration rate of the outer planet is only independent of the smoothing length if the 
sound speed is given by equation 12.71 (see next section). 

The planets are initialised on circular orbits, slowly turning on their mass in the 
first 5 orbits. We assume the re is no mass accretion , so that these remain fixed in the 
simulation (see discussion in iPeplihski et al.l l2008ai and also Sect. 12. 4p . The initial 
surface density profile is constant, and tests indicate that varying the initial surface 
density profile does not change the outcome very much. The total disc mass in the 
simulations listed in table [XT] ranges from 0.014 to 0.055 solar masses. Non-refiecting 
boundary conditions have been used throughout this chapter. 



2.3.2 Equation of state 



The outer planet is likely to undergo rapid type III mi gration, and the co-orb ital region 
will be very asymmetric. We find, in accordance with Peplihski et al.l (2008a), that the 
standard softening descript ion does not lead to convergent results (for a comparative 
study see lCrida et al.ll2009l ). Because of the massive disc, a high density spike develops 
near the planet. Any small asymmetry will then generate a large t orque, leading to 
erratic results. We follow the prescription of IPeplihski et al.l (l2008al ) and increase the 
sound speed near the outer planet since the locally isothermal model breaks down in 
the circum-planetary disc. The new sound speed is given by 



(2.7) 



where and Vp are the distance to the star and the outer planet, respectively. Qs 
and Qp are the Keplerian angular velocity and hg and hp are the aspect ratio, both of 
the circumstellar and circum-planetary disc, respectively. The parameter n is chosen 
to be 3.5, and the aspect ratio of the circum-planetary disc is hp = 0.4. The sound 
speed has not been changed in the vicinity of the inner, less massive planet as it will 
not undergo type III migration and the density peak near the planet is much lower. 

In radiative simulations (R), we go beyond the locally isothermal approximation and 
include the full thermal energy equation, which takes the generated viscous heat and 
radiative transport into account. In this radiative formulation, the sound speed is a 
direct outcome of the simulation, so we do not use equation 12.71 ( iRein et al.ll2010l ). 
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Figure 2.4: The semi-major axes (top), period ratio P2/P1 (middle), and eccentricities 
(bottom) of the two planets plotted as a function of time in dimensionless units for 
run F5 with a disc aspect ratio oi h = 0.07. In the bottom panel, the upper curve 
corresponds to the inner planet. 
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Figure 2.5: Same plot as above for simulation F4 with a disc aspect ratio of /i = 0.04. 
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2.3.3 Simulation results 



We find several possible simulation outcomes. These are indicated in the final column 
of table I2.H which gives the resonance obtained in each case with a final letter D 
denoting that the migration became ultimately divergent, resulting in a loss of the 
commensurability. We find that convergent migration could lead to a 2:1 resonance 
that was set up in the initial stages of the simulation (F3). Cases Fl ,F2,F4,F5,R1 
provide more rapid and consistent convergent migration than F3 and can attain a 3:2 
resonance directly. 

Thus avoidance of the attainment of sustained 2:1 commensurability and the effec- 
tive attainment of 3:2 commensurability required a rapid convergent migration, as pre- 
dicted by the N-body simulations and the analytic estimate (see section [2.1.2p . That is 
apparently helped initially by a rapid inward migration phase of the outer planet which 
shows evidence of type III migration. The outer planet went through that phase in all 
simulations with a surface density higher than 0.00025, and a 3:2 commensurability was 
obtained. Th us, a surface density comparable to the minimum solar nebula (MMSN, 



Havashilll98ll ). as used in simulation F3, is not sufficient to allow the planets undergo 



rapid enough type III migration. As soon as the planets approach the 3:2 commensu- 
rability, type III migration stops due to the interaction with the inner planet and the 
outer planet starts to migrate in a standard type II regime. 

This imposes another constraint on the long-term sustainability of the resonance. 
The inner planet remains embedded in the disc and thus potentially undergoes a fairly 
rapid inward type I migration. If the type II migration rate of the outer planet is slower 
than the type I migration rate of the inner planet, then the planets diverge and the 
resonance is not sustained. A precise estimate of the migration rates is impossible in 
late stages, as the planets interact strongly with the density structure imposed on the 
disc by each other. 

Accordingly, outfiow boundary conditions at the inner boundary that prevent the 
build up of an inner disc are more favourable to the maintenance of a 3:2 commensu- 
rability because the type I migration rate scales linearly with the disc surface density. 
However, those are not presented here, as the effect is weak and we stop the simulation 
before a large inner disc can build up near the boundary. 

The migration rate for the in ner planet depends on the aspect ratio h. It is decreased 
for an increased disc thickness (ITanaka et al.ll2002l ). That explains why models with a 
large disc thickness (F5,R1) tend to stay in resonance, whereas models with a smaller 
disc thickness (F2,F4) tend towards divergent migration at late times. The larger 
thickness is cons istent with radiative runs (Rl) which give a thickness of /i ~ 0.075 (see 



Rein et al.l 120101 ). 



As an illustration of the evolution of typical configuration (F5) that forms and main- 
tains a 3:2 commensurability during which the orbital radii contract by a factor of at 
least ~ 2 we plot the evolutions of the semi-major axes, the period ratio, and eccen- 
tricities in figure 12.41 We also provide surface density contour plots in figure 12.61 after 
20, 40, 60, 80, 100, 120, 140, 160 and 180 initial inner planet orbits. The eccentricity 
peak in figure 12.41 at t ~ 300 comes from passing through the 2:1 commensurability. 
At t ^ 800, the 3:2 commensurability is reached and maintained until the end of the 
simulation. The surf ace density in t his simulation is approximately 5 times higher than 
the MMSN at 1 AU flHavashilll98"lh . 

We also present an illustration of the evolution of the semi major axes and eccen- 



22 



tricities, as well as surface density plots from a simulation (F4) that does form a 3:2 
commensurability, but loses it because the inner planet is migrating too fast in figures 
12.51 and 12.71 One can see that a massive inner disc has piled up. This and the small 
aspect ratio of h = 0.04 make the inner planet go faster than the outer planet, which 
has opened a clear gap. The commensurability is lost at t ~ 1200. 

In the radiative run Rl, the initial type III migration rate is slower when compared 
to run F5 because the surface density is lower. After the the 2:1 resonance is passed, 
the outer planet migrates in a type II regime, so that the capture into the 3:2 resonance 
appears later. However, the orbital parameters measured at the end of the simulations 
are very similar to any other run that we performed. This indicates that the parameter 
space that is populated by this kind of planet-disc simulation is very generic. 

We tested the effect of disc dispersal at the late stages in our models to evolve the 
system self-consistently to the present day. In model F5 after t = 2000 we allow the 
disc mass to exponentially decay on a ti mescale of Tdis ^ 2000. This timescale is shorter 
than the photo-evaporation timescale ( lAlexander et al.l 120061). However, this scen ario 
is expected to give a stronger effect than a long timescale (ISandor fc Kleyl 120061 ). In 
agreement with those authors, we found that the dynamical state of the system does 
not change for the above parameters. At a late stage, the resonance is well established 
and the planets undergo a slow inward migration. Strong effects are only expected if 
the disc dispersal happens during the short period of rapid type HI migration, which 
is very unlikely. We observed that the eccentricities show a trend toward decreasing 
and the libration amplitudes tend toward slightly increasing during the dispersal phase. 
However, these changes are not different than what has been observed in runs without 
a disappearing disc. 

Even though we use a high surface density in our simulations, the Toomre Q pa- 
rameter is still larger than unity at the outer boundary (typically 1.8). This gives us 
confidence that we do not need to include the effects of self-gravity in the calculation. 
Furthermore, it is not expected that self-gravity plays an important role for the mi- 
gration rate of the outer planet which undergoes type III migration (M.K. Lin, private 
communication) . 
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Figure 2.6: Surface-density contour plots for simulation F5 after 20, 40, 60, 80, 100, 
120, 140, 160 and 180 inner planet orbits. At the end of the type III migration phase, 
the outer planet establishes a definite gap, while the inner planet remains embedded at 
the edge of the outer planet's gap. 
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Figure 2.7: Surface-density contour plots for simulation F4 after 20, 40, 60, 80, 100, 
120, 140, 160 and 180 inner planet orbits. The extended ID grid is not shown. 
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2.4 Other scenarios for the origin of HD45364 



In the above discussion, we have considered the situation when the planets attain their 
final masses while having a wider separation than required for a 2:1 commensurability, 
and found that convergent migration scenarios can be found that bring them into the 
observed 3:2 commensurability by disc planet interactions. However, it is possible that 
they could be brought to their current configuration in a number of different ways as 
considered below. It is important to note that, because the final commensurable state 
results from disc planet interactions, it should have similar properties to those described 
above when making comparisons with observations. 

It is possible that the solid cores of both planets approach each other more closely 
than the 2:1 commensurability before entering the rapid gas accretion phase and attain- 
ing their final masses prior to entering the 3:2 commensurability. Although it is difficult 
to rule out such possibilities entirely, we note that the cores would be expected to be 
in the super earth mass range, where in general closer comr nensurabilities than 2:1 and 
even 3:2 are found for typical type I migration rates (e.g. iPapaloizou fc Szuszkiewicz 



20051 : ICresswell fc Nelsonll2008l ). One may also envisage the possibility that the solid 



cores grew in situ in a 3:2 commensurability, but this would have to survive expected 
strongly varying migration rates as a result of disc planet interactions as the planets 
grew in mass. 

Another issue is whether the embedded inner planet is in a rapid accretion phase. 
The onset of the rapid accretio n phase (also called phase 3) occurs when the core and 
envelope mass are about equal f jPollack et al.lll996b| ). The total planet mass depends at 
this stage on the boundary conditions of the circum-planetary disc. When these allow 
the planet to have a significant convective envelope, th e transition to rapid accretion 
may not occur until the planet mass exceeds G OM^rp (jWuchterll 119931) . which is the 



mass of the inner planet (se e also model J3 of iPoUack et al.l Il996bl . and models of 



Papaloizou fc TerquemI Il999l ). Because of the above results, it is reasonable that the 



inner planet is not in a rapid accretion phase. 

Finally we remark that proto-planetary discs are beli eved to maintain tu r bulen ce in 
some part of their structure. Using the prescription of iRein fc Papaloizoul (120091 ) . we 
can simulate the turbulent behaviour of the disc by adding stochastic forces to an N- 
body simulation (see also chapter [3]). These forces will ultimately eject the planets from 
the 2:1 resonance should that form. Provided they are strong enough, this can happen 
within the lifetime of the disc, thus making a subsequent capture into the observed 3:2 
resonance possible. 

We confirmed numericall y that such cases can occu r for moderately large diffusion 
coefficients (as estimated by iRein fc Papaloizoul 120091 ) . However, this outcome seems 
to be the exception rather than the rule. Should the 2:1 resonance be broken, a planet- 
planet scattering event appears to be more likely. In all the simulations we performed, 
we find that only a small fraction of systems (1% - 5%) eventually end up in a 3:2 
resonance. 



2.5 Comparison with observations 



Table [2^2] lists the orbital parameters that ICorreia et al. ( 2008 ) obtained from their best 
statistical fit of two Keplerian orbits to the radial velocity data. The radial velocity 
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data collected so far is insufficient for detecting any interactions between the planets. 
This is borne out by the fact that the authors did not obtain any improvement in terms 
of minimising when a 3-body Newtonian fit rather than a Keplerian fit was carried 
out. However, a stability analysis supports the viability of the determined parameters, 
as these lie inside a stable region of orbital parameter space. The best-ffi solution shows 
a 3:2 mean motion commensurability, although no planet-planet interactions have been 
observed from which this could be inferred directly. It is important to note that the 
large minimum associated with the best-fit that could be obtained indicates large 
uncertainties that are not accounted for by the magnitudes of the errors quoted in their 
paper. 

There are two main reasons why we believe that the effective errors associated with 
the observations are large. First, many data points are clustered and clearly do not 
provide a random time sampling. This could result in correlations that effectively reduce 
the total number of independent measurements. Second, the high minimum value 
that is three times the quoted observational error indicates that there are additional 
effects (e.g. sunspots, additional planets, etc.) that produce a jitter in the central star 
(M. Mayor, private communication), which then enhances the effective observational 
error. 

No matter what process generates the additional noise, we can assume in a ffist 
approximation that it follows a Gaussian distribution. We then have to conclude that 
the effective error of each observation is a factor ~ 3 larger than reported in order to 
account for the minimum x^, given the quoted number of ~ 60 independent observa- 
tions. Under these circumstances, we would then conclude that any fit with \/~^ in 
the range 2 — 4 would be an equally valid possibility. 

In the following, we show that indeed a variety of orbital solutions match the ob- 
served RV data with no statisti cally significant difference when compared to the quoted 
best fit by ICorreia et al.l ( 120081 ). As one illustration, we use the simulation F5 obtained 
in section 12.31 We take the orbital parameters at a time when the orbital period of the 
outer planet is closest to the observed value and integrate them for several orbits with 
our N-body code. The solution is stable for at least one million years. There are only 
two free parameters available to fit the reflex motion of the central star to the observed 
radial velocity: the origin of time (epoch) and the angle between the line of sight and 
the pericentre of the planets. 

We can safely assume that the planet masses and periods are measured with high 
accuracy and that only the shape of the orbit contains large errors. Our best fit results 
in an unreduced a/x^ value of 3.51. According to the above discussion, this solution is 
statistically indistinguishable from a solution with ~ 2.8. 

It is possible to reduce the value of x^ even further, when assuming that the planet 
masses are not fixed. In that case, we have found a fit corresponding to \f)^ = 2.76 
where we have reduced the radial velocity amplitude by 8% (effectively adding one free 
parameter to the fit). This could be explained by either a heavier star, less massive 
planets, or a less inclined orbit. The orientation of the orbital plane has been kept fixed 
{i = 90°) in all fits. 

The results are shown in figure 12.81 The blue curve corresponds to the outcome of 
simulation F5 and we list the o rbital parameters in table 12.21 We also plot the best 
fit given by ICorreia et al.l (120081 ) for comparison (green curve). In accordance with the 
discussion given above, it is very difficult to see any differences in the quality of these 
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Correia et al. f2008) 


Simulation F5 


Parameter 


Unit 


b c 


b c 


M sin i 


L J up J 


0.1872 0.6579 


0.1872 0.6579 




[Mq] 


0.82 


0.82 


a 


m 


0.6813 0.8972 


0.6804 0.8994 


e 




0.17 ±0.02 0.097 ±0.012 


0.036 0.017 


X 


[deg] 


105.8 ± 1.4 269.5 ±0.6 


352.5 153.9 




[deg] 


162.6 ±6.3 7.4 ±4.3 


87.9 292.2 






2.79 


2.76* (3.51) 


Date 


[JD] 


2453500 


2453500 



Table 2.2: Orbital parameters of HD45364b and HD45364c from different fits. 

"Note that 137 is not well constrained for nearly circular orbits. 
''Radial velocity amplitude was scaled down by 8%. 
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Figure 2.8: Comparison of different orbital solutions. See text for a description of the 
different models. Radial velocity mea surements (red point s) and the published orbital 
solution (green curve) are taken from lCorreia et al.l (120081 ). 



fits, which indeed suggests that the models are statistically indistinguishable. 

However, there is an important diffe rence between the fits obtained from our simu- 
lations and that of Correia et al.l (2008). Our models consistently predict lower values 
for both eccentricities ei and 62 (see table [52]) • Furthermore, the ratio of eccentrici- 
ties 61/ 62 is higher than the previously reported value of 1.73. The eccentricities are 
oscillating and the ratio is on average ~ 3. 

This in turn results in a d ifferent libration pattern: the slow libration mode (see 
figure 3a in ICorreia et al.ll2008l ) that is associated with oscillations of the angle between 
the two apsidal lines is absent, as shown in figure 127101 Thus there is a marked difference 
in the form the interaction between the two planets takes place. 

It is hoped that future observations will be able to resolve this issue. The evolution 
of the difference in radial veloc ity over the n e xt few years between our fit from simu- 
lation F5 and the fit found by I Correia et al.l ( l2008l ) is plotted in figure 12.91 There is 
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Figure 2.9: Differe nce in radial velocity between the previously reported fit by 



Correia et al.l ( 120081 ) and the new fit obtained from simulation F5. 



approximately one window each year that will allow further observations to distinguish 
easily between the two models. The large difference at these dates is not caused by the 
secular evolution of the system, but is simply a consequence of smaller eccentricities. 

All runs that yield a final configuration with a 3:2 MMR are very similar in their 
dynamical behaviour, as found for model F5. They are characterised by an antisym- 
metric {vu2 — OTi ~ vr) state with a relatively small libration of 0i and 02. This is a 
generic outcome of all convergent disc-planet migration scenarios that we performed. 



2.6 Conclusions 

The planets in the multi-planetary system HD45364 are most likely in a 3:2 mean motion 
resonance. This poses interesting questions on its formation history. Assuming that 
the planets form far apart and migrate with a moderate migration rate, as predicted 
by standard planet formation and migration theories, the most likely outcome is a 2:1 
mean motion resonance, contrary to the observation of a 3:2 MMR. 

In this chapter, we investigated a possible way around this problem by letting the 
outer planet undergo a rapid inward type III migration. We presented an analytical es- 
timate of the required migration rate and performed both N-body and hydrodynamical 
simulations. 

We find that it is indeed possible to form a 3:2 MMR and avoid the 2:1 resonance, 
thus resembling the observed planetary system using reasonable disc parameters. Hy- 
drodynamical simulations suggest that the system is more likely to sustain the resonance 
for high aspect ratios, as the migration of the inner planet is slowed down, thus avoiding 
divergent migration. 

Finally, we used the orbital configuration found in the hydrodynamical formation 
scenario to calculate a radial velocity curve. This curve was then compared to obser- 
vations and the resulting fit has an identical value to the previously reported best 
fit. 
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Figure 2.10: Evolution of the resonant angles 0i = 2Ai — 3A2 + 02 = 2Ai — 3A2 + ^^2, 
zui — W2 for simulation F5 using the parameters given in table 12.21 The angles are 
measured in radians. 

Our solution is stable for at least a million years. It is in a dynamically different 
state, both planets having lower eccentricities and a different libration pattern. This is 
the first time that planet migration theory can predict a precise orbital configuration 
of a multiplanetary system. This might also be the first direct evidence for type III 
migration if this scenario turns out to be true. 

The system HD45364 remains an interesting object for observers, as the differences 
between the two orbital solutions can be measured in radial velocity within a couple of 
years. 
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Chapter 



The survival of mean motion 
resonances in a turbulent disc 



Every body continues in its state of rest, or of 
uniform motion in a right line, unless it is compelled 
to change that state by forces impressed upon it. 

Isaac Newton, Principia Mathematica, 1687 



As shown for the system HD45364 in chapter [21 resonant configurations can be estab- 
hshed as a result of dissipative forces acting on the planets which lead to convergent 
migration. The previous chapter assumes a laminar disc. However, in general the disc is 
thought to be turbulent because, as discussed in more detail in chapter [H the molecular 
viscosity is not large enough to drive the observed outward angular momentum flux. 
The most promising mechanism t o drive turbulence is the magneto- rotational instabil- 
ity, MRI (IBalbus fc Hawlevlll99l[). although other possibilities are still being discussed 
(see e.g. iLesur fc Papaloizoull2009[ ). 

In this chapter, we study the effects of density perturbations in the disc resulting 
from any kind of turbulence on planetary systems. We do not solve the full 3D magneto 
hydrodynamics equations, but rather assume a parametrisation for the forces acting 
on planets. This has two main advantages. First, we can understand the relevant 
dynamical processes much more easily. Second, there is a large diversity of simulations 
of the MRI in the literature which exhibit different diffusion coefficients, varying by 
several orders of magnitude. By parametrising the forces, we avoid these difficulties 
altogether. 

Using this approach, we study the response of both a single planet and planetary 
systems in a 2:1 mean motion commensurability to stochastic forcing. The forcing could 
for example result from MRI turbulence. We first develop an analytic model from first 
principles. We can isolate two distinct libration modes for the resonant angles which 
react differently to stochastic forcing. 

Systems are quickly destabilised if the magnitude of the stochastic forcing is large. 
The growth of libration amplitudes is parametrised as a function of the diffusion co- 
efficient and other relevant physical parameters. We also perform numerical N-body 
simulations with additional stochastic forcing terms to represent the effects of disc tur- 
bulence. These are in excellent agreement with the analytic model. 
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Stochastic forcing due to disc turbulence may have played a role in shaping the 
configurations of observed systems in mean motion resonance. It naturally provides a 
mechanism for accounting for the HD128311 system, for which the fast mode librates 
and the slow mode is apparently near the borderline between libration and circulation. 



3.1 Basic equations 



Papaloizoul l2003[ ) 



ie planet moving 


in a 


Snellerove et al. 


2001; 



E . -JJi ,3.1, 

dH dH\ 



3X ' Or.) P-^' 

Here the angular momentum of the planet is L and the energy is E. For a Keplerian 
orbital motion around a central point mass M we have 



m 



^GMa{l - e2) and (3.5) 

GMm , , 

E = . (3.6) 

where G is the gravitational constant, a the semi-major axis and e the eccentricity (see 
appendix For an elliptical orbit around a point mass, equation 13.31 results in a linear 
growth of the mean longitude 

\ = n{t- to) + U7, (3.7) 

where n is the mean motion, to being the time of periastron passage and w being the 
longitude of periastron. All other equations of motion are trivial as the right hand side 
equates to zero. 



3.1.1 Additional forcing of a single planet 

In order to study phenomena such as stochastic forcing, we need to consider the effects of 
an additional external force per unit mass F, which may or may not be described using 
a Hamiltonian formalism. However, as may be seen by considering general coordinate 
transformations starting from a Cartesian representation, the equations of motion are 
linear in the components of F. Because of this we may determine them by considering 
forces of the form F = [F^, Fy) for which the Cartesian components are constant. 
Having done this we may then suppose that these vary with coordinates and time in an 
arbitrary manner. Following this procedure we note that when F, as in the above form 
is constant, we can derive the equations of motion by replacing the original Hamiltonian 
with a new Hamiltonian defined through 

H ^ H - m{F^ X + Fy y) = H - m{Y ■¥) . (3.8) 
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The additional terms proportion al to the components of F c orrespond to the Gaussian 
form of the equations of motion ( iBrouwer fc Clemencdll96lh . 

The various derivatives involving r can be calculated by elementary means and 
expressed in terms of E, L, A and w. One thus finds additional contributions to the 
equations of motion 13.11 - 13. 4[ indicated with a subscript F, in the form 



Lf 
Ep 

Wp 



d_ d_ 

^dX dw 

d , 
mn— (r ■ F) 

OA 

-m- (r . F) 



na e 



—m 



(r ■ F) = m (r X F) 

m (v ■ F) 

or equivalent ly 
1 



1 

r F 



a 



sin f — Fr cos / 



ZO p 



2an 
GM 



(r-F) 



(3.9) 

(3.10) 
(3.11) 

(3.12) 
(3.13) 



where the true anomaly / is defined as the difference between the true longitude and 
the longitude of periastron, f = 6 — w (see also appendix . Note that from equation 
I3.1UI we obtain 

2an 2{Fresmf + Fe{l + ecosf)) 



ap = 

Sn nV 1 - 

and from equation 13.91 together with equation 13.101 we obtain 



ep 



L{2EL + LE) 



(3.14) 



(3.15) 



(3.16) 



In the limit e ^ 1 this becomes (ignoring terms 0(e) and smaller) 

ep = Fr — sin/ + Fe — 2 cos/. 

an an 

Furthermore in this limit we may replace f hj f = X — to = n(t — to) . 

Note that that the above formalism results in equation l3.12l for lijp, which diverges 
for small e as 1/e. This comes from the choice of coordinates used and is not associated 
with any actual singularity or instability in the system. This is readily seen if one 
uses h = e since, and k = ecoscc as dynamical variables rather than e and zu. The 
former set behave like Cartesian coordinates, while the latter set are the corresponding 
cylindrical polar coordinates. When the former set is used, potentially divergent terms 
oc 1/e do not appear. This can be seen from equations 13.121 and 13.161 which give in the 
small e limit 



hp = —Fr — sinA + Fe — 2 cos A 

an an 

kp = Fr — sinA + Fe — 2 cos A. 

an an 



(3.17) 
(3.18) 



Abrupt changes to zu may occur when h and k pass through the origin in the [h, k) 
plane. But this is clearly just a coordinate singularity rather than a problem with the 
physical system which changes smoothly on transition through the origin. The abrupt 
changes to the zu coordinate occur because very small perturbations to very nearly 
circular orbits produce large changes to this angle. 
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3.1.2 Multiple planets 



Up to now we have considered a single planet. We now generalise the forma lism so that 
i t app lies t o a system of two planets. We fo llow closely the discussions in IPapaloizou 



( 2003 ) and Papaloizou &: Szuszkiewic3 fl2005 ). Excluding stochastic forcing for the time 
being, we start from the Hamiltonian fo rmalism describing their mutual interactions 
using Jacobi coordinates fISinclaijflQTsI l In this formalism the motion of the outer 
planet is described around the centre of mass of the star plus the inner planet rather 
then the star alone. Thus, the radius vector r2, of the inner planet of mass m2 is 
measured from the star with mass M and that of the outer planet, ri, of mass mi is 
referred to the centre of mass of M and m2. The reduced masses of m2 and mi are 
then defined as fi2 = ■ M/ (M + m2) and /xi = mi ■ (M + m2)/(M + m2 + mi). We 
also define M2 = M + m2 and Mi = M + m2 + mi. From now on we consistently adopt 
subscripts 1 and 2 for coordinates related to the outer and inner planets respectively. 
Note that this is different from the previou s chapter 



The required Hamiltonian, is given by (IMurray fc Dermottll2000l p.440ff): 



H 



1 2^1 |. |2 GM2fl2 



GMlTLi 017111712 



I roil 



|ri2| 



(3.19) 



Here is the relative position vector between the objects with subscript i and j, where 
denotes the star. The Hamiltonian can be expressed in terms of Ei,Li,Wi,Xi,i = 
1,2 and the time t. The energies are given by Ei = —GifiiMi/{2ai) and the angular 
momenta Li = ^i^JGMiai{l — ef) with and Cj denoting the semi-major axes and 
eccentricities respectively. The mean motions are rii = we 
replace Mj by M and /Zj by m^ from now on. 

The interaction Hamiltonian (the term oc mim2 in equation I3.19P can be expanded 
in a Fourier series involving li near combinations of the t hree angular differences A, — 



zji, i = 1,2 and zui — zu2 (e.g. iBrouwer fc ClemencdllQGll ). Near a first order p + 1 : p 



{p + l)Ai - p\2 - ^2-, and 



resonance, we expect that both 0i 
jpA^ — w^ , will be slowly varying. Following standard practice (see e.g. lPapaloizoull2003 



Papaloizou fc Szuszkiewiczl 120051 ) only terms in the Fourier expansion involving linear 
combinations of 0i and 02 as argument are retained because only these are expected to 
lead to large long term perturbations. 

The resulting Hamiltonian may then be written in the general form H = Ei + E2 + 
H12, where the interaction Hamiltonian is given by 



H 



12 



Gmim2 
ai 



(-, 



Ci, 62 ) cos(/0i + A;0 



2J, 



(3.20) 



where in the above and similar summations below, the sum ranges over all positive and 
negative integers (fc, /) and the dimensionless coefficients C^^i depend on 61,62 and the 
ratio of semi-major axes ai/a2 only. 



Equations of motion for two planets 



The equations of motion for each planet can now be easil y derived, that take into ac- 



count the contributions due to their mutual interactions (see Papaloizoull2003 : Papaloizou fc Szuszkiewi' 



20051 ) and contributions from 13.91 - 13.131 The latter terms arising from external forcing 
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are indicated with a subscript F. We obtain to lowest order in the perturbing masses: 
drii 



dt 

dn2 
~dt 



3(p +1^^^ ^ C,,(^ + Ml^. + ^fe) + (^) ^ (3.21) 



3pnlmia2 
Ma'i 



dn2 
~dt 



(3.22) 



dci 
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_l_ / dei 

\dt y p 



aie2M 
l+p{k + l) [ 1 



Ck,i sin(/0i + A;02) 
rfe2 



1-ei 



+ 



dt 



(3.23) 



(3.24) 



Here 



and 



: (p + l)ni - pn2 - 5^(i^fe,i + cos(/</.i + k<P2) + (^^) (3.25) 
= (p+l)ni-pn2-5^(i^fe,/ + Ffc,;)cos(/0i + A;02)+ (^^^ • (3.26) 
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(3.27) 



(3.28) 



(3.29) 



Note that 02 — 0i = "^2 — t^i = C is the angle between the two apsidal lines of the two 
planets. We also comment that, up to now, we have not assumed that the eccentricities 
are small and that, in addition to stochastic contributions, the external forcing terms 
may in general contain contributions from very slowly varying disc tides. 
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Modes of libration 



We first consider two planets in resonance with no external forces, to identify the 
possible libration modes. In the absence of external forces, equations 13.211 - 13.261 can 
have a solution for which Oj and Cj are constants and the angles 0i and 02 are zero. 
In general, other values for the angles may be possible but such cases do not occur for 
the numerical examples presented below. When the an gles are zero, equ ations 13.251 and 



13.261 provide a relationship between ei and 62 (see e.g. iPapaloizod l2003l ) . 

We go on to consider small amplitude oscillations or librations of the angles about 
their equilibrium state. Because two planets are involved there are two modes of os- 
cillation, which we find convenient to separate and describe as fast and slow modes. 
Assuming the planets have comparable masses, the fast mode has a libration frequency 
oc yjm and the slow mode has a libration frequency oc m. These modes clearly separate 
as the planet masses are changed while maintaining fixed eccentricities. 

Fast mode 

To obtain the fast mode we linearise equations 13.211 - 13.261 and neglect second order 
terms in the planet masses. This is equivalent to neglecting the variation of Dk i,Ek^i 
and Fk i in equations 13.251 and 13.261 which then require that (pi = 02 very nearly for this 
mode. Noting that for linear modes of the type considered here and in the next section, 
equations 13.211 - 13.241 imply that and Cj are proportional to linear combinations of 
the librating angles, differentiation of either of equations 13.251 or 13.261 with respect to 
time then gives for small amplitude oscillations 

+ = 0, (1 = 1,2), (3.30) 



where 

2 3p^n2m2 



(l + ^).J2CUk + l?, (3.31) 



Jif — 

M 

and we have used the resonance condition that (p + l)ni = pn2 which is satisfied to 
within a correction of order ^JmxjM. Note that for this mode the fact that 0i = 02 
very nearly, implies that Wi — W\ is small. Thus that quantity does not participate in 
the oscillation. 

Slow mode 

Now, we look for low frequency librations with frequency oc m\. Equations 13.251 and 
13.261 implv that, for such oscillations, to within a small relative error of order uii/M, 
{p + l)ni = pn2 throughout. 

Equations 13.211 and 13.221 then imply that the two angles are related by 02 = /30i, 
where (3 = —J2^k,i{k + l)l)/{J2Ck,i{k + l)k). Subtracting equation 13.261 from equation 
I3.25[ differentiating with respect to time and using this condition results in an equation 

for C = 02 - 01 = -072 - lUi 

^ + ^IC = 0, (3.32) 

where 



2_ nim2^1 - ej dW nzazmiy/l - e| dW 
"^'^ - e,M{l-/3) 9^7 + aie2M(l-/3) ^^'^^^ 
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with 



Q!2 = ^ Ck,l 



k-{p+l){k + l) [l-Jl-ej 



{k(5 + l), 



l + p{k + l)[l- a/1 -e: 



{kl3 + /) 



and 



W 



(nia2m2\/l - e| d n2rax\/\ - ef d \ sr^ ^ 
aiezM ^^M 9^ J ' ^ 



Although the expressions for the mode frequencies are comphcated, the fast frequency 
scales as the square root of the planet mass and the slow frequency scales as the planet 
mass. Both frequencies scale with the characteristic mean motion of the system. 

Librations with external forcing 

When external forcing is included source terms appear on the right hand sides of equa- 
tions [230] and 12321 We assume that the forcing terms are small so that terms involving 
products of these and both the libration amplitudes and the planet masses may be 
neglected. Then, in the case of the slow mode, repeating the derivation given above 
including the forcing terms, we find that equation 13.321 becomes 

^ + ^f.C = ^^{w2F - -^xf) ■ (3.34) 

The quantities Wip are readily obtained for each planet from equation I3.12[ From 
this we see that for small eccentricities, Wip oc l/cj, indicating large effects when Cj 
is small. As already discussed, this feature arises from a coordinate singularity rather 
than physically significant changes to the system. 

A similar description may be found for the fast mode. In this case, neglecting terms 
of the order of the square of the planet masses or higher, one may use equations 13.251 
and 13.261 to obtain an equation for 0i in the form 

+ = — \^xF^ + (P + y)nxp - pn2F. (3.35) 

Equations 13.351 and 13.341 form a pair of equations for the stochastically forced fast 
and slow modes respectively. This mode separation is not precise. However, it can be 
made so by choosing appropriate linear combinations of the above modes. Numerical 
results confirm that (pi predominantly manifests the fast mode and ( the slow mode, 
so we do not expect such a change of basis to significantly affect conclusions. 

We further comment that because (pip contains W2p but not wip, for small ec- 
centricities there are only potential forcing terms oc l/e2 that occur when forcing is 
applied to the inner planet. As this planet has the larger eccentricity for the situations 
we consider, small eccentricities are not found to play any significant role in this case. 

Each mode responds as a forced oscillator. The forcing contains a stochastic compo- 
nent which tends to excite the respective oscillator and ultimately convert libration into 
circulation. But we stress that the above formulation as well as developments below 
assume small librations, so we may only assess the initial growth of oscillation ampli- 
tude. However, inferences based on the structure of the non linear governing equations 
and an extrapolation of the linear results enable successful comparison with numerical 
results. 
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3.1.3 Stochastic forces 



We assume that turbulence causes the external force per unit mass {Fr, Fq) acting on 
each planet to be stochastic. For simplicity we shall adopt the simplest possible model. 
Each component of the force satisfies the relation 

F.mit') = {F^)g{\t-t'\) (3.36) 

where the auto-correlation function g[x) is such that g{x)dx = Tc, where Tc is the 
correlation time and ^^/{Ff) is the root mean square value of the i component. Note 
that an ensemble average is implied on the left hand side. Again, for simplicity we 
assume that different components acting on the same planet as well as components 
acting on different planets are uncorrelated. The root mean square values as well as 
may in general depend on t, but we shall not take this into account here and simply 
assume that these quantities are constant and the same for each force component. 

The stochastic forces make quantities they act on undergo a random walk. Thus, if 
for example A = Fi for some quantity^ A, the square of the change of A occurring after 
a time interval t is given by 

"t ft 



(A A) 2 =11 Fi{t')F,{t")dt'dt" 
Jo Jo 

= [ [ {Fl)g{\t' -t"\)dt'dt" ^ Dit. (3.37) 
Jo Jo 

Here we take the limit where t/r^ is very large corresponding to an integration time of 
many correlation times and Z)j = 2 {Ff) is the diffusion coefficient. 

When the evolution of a stochastically forced planetary orbit is considered, it is 
more appropriate to consider a model governing equation for A of the generic form 

A = FiSm{nt), (3.38) 

where we recall that 27r/?T, is the orbital period (but note that a different value could 
equally well be considered). We note in passing that, by shifting the origin of time, 
an arbitrary phase may be added to the argument without changing the results given 
below. One readily finds that equation 13.371 is replaced by 

{AAy = f f Fi{t')Fi{t")sm{nt')sm{nt")dt'dt" 
Jo Jo 
t pt 

(F^-) g{\t' - ("I) sin(n«') sHnt")dt'dt" 



- (3.39) 

where we introduce a new correction factor 7(n) = g{x) co?>{nx)dx. If riTc ^ 1, 
corresponding to the correlation time being much less than the orbital period, then 
7 — > 1. For larger Tc, 7 < 1 gives a reduction factor for the amount of stochastic 



^Constants or slowly varying quantities originally multiplying Fi may be absorbed by a redefinition 
of A and so do not affect the discussion given below. 
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diffusion. For example if we adopt, an exponential form for the auto-correlation function 
such that 



g{\t'-t"\)=exp 



t' - t" 



we find 

^ = 1 \ 2 ' (3-40) 

and for the purposes of simplicity and comparison with numerical work we shall use 
this from now on. 

Stochastic forcing of an isolated planet 

In the case of a stochastically forced single planet we may obtain a statistical estimate 
for the characteristic growth of the orbital parameters as a function of time by directly 
integrating equations EUHl EII21 13. 16^ l3.l71 and [XTS] with respect to time. We may then 
apply the formalism leading to the results expressed in generic form through equations 
13. 371 - 13^^ to obtain estimates for the stochastic diffusion of the orbital elements in the 
limit of small eccentricity in the form 

(Aa)2 : 
(Aef : 



4 


(3.41) 


Dt 


(3.42) 




(3.43) 


2.5 -fDt 


(3.44) 




(3.45) 




(3.46) 



As mentioned before, the 1/e^ dependence of (Acc)^ which arises from the coordinate 
singularity does not appear for (A/i)^ and (Afc)^. Note that the definitions of h and k 
imply, consistently with the above, that 

(A/i)2 = {Aef sin^ zu + e^{ Aw fcos'^vj (3.47) 
(AA;)^ = (Aefcos^w + e'^iAwfsm^w. (3.48) 

The evolution of the longitude 

Of special interest is also the evolution of the planet's mean longitude. In observations 
of Saturn's rings, for example, it is much easier to measure a shift of the longitude of an 
embedded moonlet compared to a direct measurement of the change in its semi-major 
axis. The time derivative of the mean motion is given by 



3 GM. 3Fe ,^ .^^ 

n = --J-^a = -^. (3.49) 
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Assuming a nearly circular (e <^ 1) orbit, the mean longitude is thus given by the 
double integral 

X{t) = ! {n + ^n{t')) dt' (3.50) 
Jo 



= nt- / ' dfdt. (3.51) 

Jo Jo 

The root mean square of the difference compared to the unperturbed orbit is 

(AA)' = ((A(t) -nt)')°'' (3.52) 

jPe{t }J^e{t J dt""dt"'dt"dt' (3.53) 



^0 ^0 ^0 

t ft' rt rt 



a? 



^^-^ f f f f g{\t" -t""\)dt""dt"'dt"dt' (3.54) 
Jo Jo Jo Jo 



(3.55) 



^-^^e = \-,e = {Anf'l (3.56) 
2 a"' 6 

where we have taken the limit of t ^ Tc in the last line. We also assume an exponentially 
decaying auto-correlation function g because the integral in equation 13.541 can not be 
solved analytically for an arbitrary function g. Note that on short timescales, other 
terms in equation 13.561 may dominate the growth of A. 

Instead of a stochastic migration, let us also assume a laminar migration with con- 
stant migration rate Ta, given by 

Ta = - = 3.57 
a 2n 

As above, we can calculate the longitude as a function of time as the following double 
integral 

A(t) = / {n + An{t')) dt' (3.58) 



nt- / dt"dt' = nt- -—t^ 3.59 

Jo Jo 2r, ATa ^ ' 



and thus 



— -t^= (An) - 



{AXy = ^—t*={Any-. (3.60) 



a 



Note that equation 13. 60l is an exact solution, whereas equation I3.56l is a statistical quan- 
tity, describing the root mean square value in an ensemble average. In the stochastic and 
laminar case, (AA)^ grows like ~ t^ and ~ t^, respectively. This allows to discriminate 
the different mechanisms by observing AA over an extended period of time. 
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Stochastic variation of the resonant angles in the two planet case 

We now consider the effects of stochastic forcing on the resonant angles. The expressions 
are more complicated than in the previous section because there are more variables 
involved. However, we basically follow the same formalism. 

While the libration amplitude is small enough for linearisation to be reasonable, the 
evolution is described by equations 13.341 and 13.351 These may be solved by the method 
of variation of parameters. Assuming the amplitude is zero at i = 0, the solution of 
equation 13.341 is given by 

r ■ ( +\ S(;cos{uJist) St;sm{ujist) 

C = sm[uist) / — ^ at — cos{uist) / — ^ at, (3.61) 

Jo <^ls Jo ^Is 

where Sc_ = d{w2F — ^iF)/dt. Equation 13.611 mav be regarded as describing a harmonic 
oscillator whose amplitude varies in time such that the square of the amplitude after a 
time interval t is given by 

^ / s,sinM \ s,cosM \ ' g^) 



Iq ^Is / V-'O ^Is 

The corresponding expression for equation 13.351 is 

^ / S,Mu,jt) \ ^ ^ / S,cos{u,,t) \ ' ^3 
\Jo ^if J \Jo ^if J 

where S(f, = d{(f)iF)/dt + {p + l)riiF — pn2F- 

We now evaluate the expectation values of these using the formalism of the previous 
section. For simplicity we specialise to the case when stochastic forces act only on the 
outer planet. This is also what has been considered numerically later. Taking equation 
I3.62[ we perform an integration by parts, neglecting the end point contributions as 
these are sub-dominant (they increase less rapidly than t for large t), to obtain 

(A^)^ = 'd7iFsm{uist)dt^ ^ '^'^^^^'■'^^^^^ ' (3.64) 

In dealing with equation 13.631 we neglect (pip in because after integration by parts 
this leads to a contribution on the order oJij/n smaller than that derived from hip- 
Thus, we simply obtain 

(A0i)2 _ f hipsm{ujift) ^ / TiiF cosjuift) 



{p+iy \Jo ujif J Vio ^if 

We again follow the procedures outlined in the previous section to obtain the final 
growth estimates 



(AeieasinC)^ = 2.5 ^ I and (3.65) 



2ajnj 
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where the corrections factors are given by 



7/ = ^^-r-, and (3.67) 
_ 1 1 

~ l + (ni+a;z,)V2 + l + (ni-a;,,)2rr ^^'^^^ 

The correction factors 7/ and 7^ compare the relevant timescales of the system (here, 
the hbration periods) to the correlation time of the stochastic force. 



Growth of libration amplitudes 

Equations 13.651 and 13.661 express the expected growth of the resonant angle libration 
amplitudes as a function of time. These expressions can be simply related to those 
obtained for a single planet. Thus equations 13.411 and 13.661 applied to the outer planet 
imply that 

(A0i)V(Aai)2 = 9(p + l)X7//(4aH)- 

As we are interested in the case p = 1, the width of the libration zone is ~ aiUif/ni, 
we see that the time for (A0i)^ to reach unity is comparable for the semi-major axis 
to diffuse through the libration zone. 

For fixed eccentricity of the inner planet 62, equation 13.651 indicates that C, being 
the angle between the apsidal lines of the two planets, diffuses in the same way as for an 
isolated outer planet subject to stochastic forces. Thus, in the small amplitude regime, 
the way this diffusion occurs would appear to be essentially independent of whether the 
planets are in resonance. 

An important consequence of equation 13.651 is the behaviour of C for small ei. The 
latter quantity was assumed constant in the analysis. As already discussed earlier, 
abrupt changes to C, are expected when the eccentricity of one planet gets very small. 
Then, even an initially small amplitude libration is converted to circulation. Thus if 
Ci is small, then equation 13.651 indicates that a time t ~ 'ia\n\e\/ {bD'^sC'^)i is required 
to convert libration to circulation. This can be small if Ci is small. Even if ei is not 
small initially, it is important to note that it also undergoes stochastic diffusion (see 
equation I3.43P as well as oscillations through its participation in libration. Should e\ 
attain very small values through this process, then from equation 13.651 we expect the 
onset of a rapid evolution of C,. This is particularly important when e\ already starts 
from relatively small values. 

In fact, application of equations 13.651 and 13.661 to the numerical examples discussed 
below, adopting the initial orbital elements, indicate that the diffusion of C, is signifi- 
cantly smaller than 0i unless ei starts out with a very small value. This would suggest 
that 01 reaches circulation before C- However, this neglects the coupling between the 
angles that occurs once the libration amplitudes become significant. It is readily seen 
that it is not expected that (pi could circulate while C, remains librating as it was initially. 
One expects to recover standard secular dynamics for C, from the governing equations 
13.211 - l3^25| when these are averaged over an assumed (pi circulating with constant 0i. 
As a libration of the initial form would not occur under those conditions, we expect, 
and find, large excursions or increases in the libration amplitude of C, to be correlated 
with increases in the libration amplitude of (pi. This in turn increases the oscillation 
amplitude of the eccentricity ei, allowing it to approach zero. The consequent rapid 
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evolution of ( then enables it to pass to circulation. Thus, the breaking of resonance 
is ultimately found to be controlled by the excitation of large amplitude librations for 
01, which induce ( to pass to circulation somewhat before 0i itself. 



3.2 Numerical simulations 

We have performed numerical simulations of one and two planet systems that allow for 
the incorporation of additional stochastic forces with the properties described above. 
These forces provide a simple prescription for estimating the effects of stochastic grav- 
itational forces produced by density fluctuations associated with disc turbulence. The 
results have been obtained using the N-body code described in chapter [21 We have 
checked that results are converged and do not depend on the integrator used. 

We first discuss the expected scaling of the stochastic forces with the physical pa- 
rameters of the disc and their implementation in N-body integrations. In order to clarify 
the physical mechanisms involved, and to check the analytic predictions for stochas- 
tic diffusion given by equations 13.411 - 13.451 we consider simulations of a single planet 
undergoing stochastic forcing first. We then move on to consider two planet systems 
with and without stochastic forcing. We focus on the way a 2:1 commensurability is 
disrupted and highlight the various evolutionary stages a system goes through as it 
evolves from a state with a strong commensurability affecting the interaction dynam- 
ics, to one where the commensurability is completely disrupted. We find that in some 
cases a strong scattering occurs. We consider a range of different planet masses and 
eccentricities (see table 13. ll) . 



3.2.1 Stochastic forces 



In order to mimic the effects of turbulence, for example produced by the MRI, it is 
necessary to calibrate these forces with reference to MHD simulations. As described 
above, the basic parameters characterising the prescription for stochastic forcing that 
we have implemented are the root mean square value of the force components per unit 
mass (in cylindrical coordinates) \/jFf) and the auto correlation time Tc. 

From our analytic considerations, we conclude that stochastic forces make the or- 
bital parameters undergo a random walk that is dependent on the force model primarily 
through the diffusion coefficient D = 2 {Ff) Tc but also through the dimensionless cor- 
rection factors 7,7s and 7/. 

For planets under the gravitational influence of a proto-planetary disc, the natural 
scale for the force per unit mass componen ts, Fj, is Fn(r) = 7r(jE(r)/2, w here S is the 
characteristic disc surface density (see e.g. iPapaloizou fc TerquemI 120061 ). Fo(r) is the 



gravitational force per unit mass due to a small circular disc patch of radius rs at a 
distance y/2 from its centre assuming that all its mass is concentrated there. See 
also figure 13.11 This is independent of r^. The natural correlation time is the inverse 
of the orbital angul ar frequency Tr n = We adopt a minimum mass solar nebula 

model (MMSN, see I Weidenschillind 1 1 9 77bl ) with 



4200- 



cm^ 
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1 AU 



-3/2 



(3.69) 



Planet 



Figure 3.1: The planet is embedded in a proto-planetary disc with random density 
fluctuations. The force exerted by both clumps in the picture on the planet has the same 
strength. It is independent of the clump size in a two dimensional system, assuming 
that the clumps are randomly distributed. 



This provides a natural scale for D as a function of the local disc radius and the central 
stellar mass through^ 



2CF^Tc,o = 1.95 C 



cm 



1 AU 



1 M 







(3.70) 



where C is a dimensionless constant. There are several factors which contribute to the 
value of C: 



The density fluctu ations found in MRI simulations Sp/p 
the order 0.1 (e.g. iNelsonI l2005l ). 



are typically of 



The presence of a dead zone in the mid plane regions of the disc, where the MRI 
is not active, has been found to cause reductions in the r nagnitude of Fn b y one 
order of magnitude or more, as compared to active cases flOishi et al.l 120071 ). 



Massive planets open a gap in the disc. lOishi et al.l (120071 ) found that most of 
the contribution to the stochastic force comes from density fluctuations within a 
distance of one scale-height from the planet. When a gap forms, this region is 
cleared of material leading to a substantial decrease in the magnitude of turbulent 
density fluctuations. Consequently Fq should be reduced on account of a lower 
ambient surfa ce density. A factor of .1 seems reasonable, although it might be 
even smaller (jde Val-Borro et al.ll2006l ). 



The correlation time Tr is found to be slightly shorter thari an orbital period, 
namely 0.5 fi"^ fiNelson fc PapaloizoulEoo3 lOishi et al.ll2007h . 



If it is appropriate to include reduction factors to account for all of the above effects, 
one flnds C = 5 ■ 10"'' and we expect a natural scale for the diffusion coefficient to be 
specified through 

o„^ro-HA.r"(^V''''4- (3.71) 



1 AU 



IMq 



^The scaling in Rein fc Papaloizoul (|2009ll is incorrect due to a missing factor of tt. The value in 
this thesis has been corrected. 
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Planet 


m (Mj) 


a (AU) 


P (days) 


e 


c 




GTS7fi p 


0.790 


0.131 


30.46 


0.263 


10° 




b 


2.530 


0.208 


60.83 


0.031 






GTS7fi T.M r 


0.13 


0.131 


30.46 


0.263 


10° 




b 


0.42 


0.208 


60.83 


0.031 






vjrJOiU kJJ_/ 


0.013 


0.131 


30.46 


0.263 


10° 




b 


0.042 


0.208 


60.83 


0.031 






G T87fi F, p 

VJ J O 1 U 11/ 


0.0013 


0.0131 


30.46 


0.263 


10° 




b 


0.0042 


0.208 


60.83 


0.031 






76 LM HE c 


0.13 


0.131 


30.46 


0.41 


10° 




b 


0.42 


0.208 


60.83 


0.09 




GJ 


876 SE HE c 


0.013 


0.131 


30.46 


0.41 


10° 




b 


0.042 


0.208 


60.83 


0.09 






HD128311 b 


1.56 


1.109 


476.8 


0.38 


58° 




c 


3.08 


1.735 


933.1 


0.21 





Table 3.1: Parameters of the model systems cons idered. Orbital elements of the system 
GJ876 are based on a fit by iRivera et al.l (120051 ). LM, SE and E indicate that planet 
masses have scaled down by a factor of 6, 60 and 600, respectively. This is equivalent to 
a Jovian mass planet, a super-Earth and Earth mass planet around a solar mass star. 
HE indicates higher eccentricities than in the observ ed system. Orbital elements of the 
system HD128311 are taken from lVogt et al.l (120051 ) . 



The same value of D may be equivalently scaled to the orbital parameters of the planets 
without reference to the disc by writing 

= ^'^■l" (lAu) (im;) X-Wkf-)—- 

Thus, a value = 10~^ in cgs units corresponds to a ratio of the root mean square 
stochastic force component to that due to the central star of about ~ 10~^ for a central 
solar mass at 1 AU. It is a simple matter to scale to other locations. 

Of course we emphasise that the value of this quantity is very uncertain, a situa- 
tion that is exacerbated by its proportionality to the square of the magnitude of the 
stochastic force per unit mass. For this reason we perform simulations for a range of 
D, covering many orders of magnitude. 



3.2.2 Numerical implementation of correlated noise 

The procedure we implemented, uses a discrete first order Markov process to generate 
a correlated, continuous noise that is added as an additional force. This statistical 
process is defined b y two parame ters, the root mean square of the amplitude and the 
correlation time Tc ( lKasdirull995l ). It has a zero mean value and has no memory. This 
has the advantage that previous values do not need to be stored in memory. The 
auto-correlation function decays exponen tially and thus rn imics the auto-correlation 
function measured in MHD simulations by lOishi et al.l (12007! ) . It is also consistent with 
the assumptions made in the previous sections. 



45 



The procedure for generating noise in 
tau and root mean square value of 1 is as 
pseudo code is 

X *= exp( -dt / tau) 

X += Xi * sqrt( 1 - exp( 2 * dt 

where Xi is a random variable with norma 



a variable x with an auto-correlation time 
follows. The update of x from t to t+dt in 

/ tau ) ) 

1 distribution and a variance of unity. 



3.2.3 Stochastic forces acting on a single planet 

We first investigate the long-term effect of stochastic forces on a single isolated planet. 
The initial orbital parameters are taken to be the observed parameters of GJ876 b (see 
table 13. Ij) . In this simulation, we use the reduced diffusion coefficients 

D = 8.2-10-'^, and D = 8.2 ■ lO'^^, 
s"* s"* 

which can be represented by a correlation time of half the orbital period and a specific 
force with root mean square values ^/{F^ = 4.05 ■ lO'^ cm/s^, and ^/{F^ = 4.05 • 
10~^ cm/s^, respectively (see also equation I3.72|) . In that case the correction factor 7 
evaluates to 7 = 0.0866. 

The growth rates of the root mean square values for Aa, Aw and Ae in 200 simu- 
lations are plotted in figure 13.21 We also plot the analytic predictions (equations 13. 4H 
13.451 and I3.43P as solid lines. The numerical simulations and the analytic model are in 
excellent agreement and scale as expected. We have performed simulations for a variety 
of diffusion coefficients and get results that are fully consistent in all cases. 

We also plot the analytic growth rates coming from equations 13.451 and 13.431 while 
setting 7 = 1, which corresponds to uncorrelated noise. The resulting curves (dashed 
lines) do not agree with simulation results. This shows clearly that corr elation effects 



are im portant, and have to be taken into account. The model presented by I Adams et al.l 



(I2OO8I ) makes use of such an uncorrelated force (random kicks), leading to an overesti- 
mated diffusion coefficient by almost a factor ~ 16. Note that only orbital parameters 
involving the geometry of the orbit are affected, i.e., not the semi- major axis. 

In figure 13. 3[ we plot the root mean square of the difference in longitude. One can 
clearly seen that (AA)^ does grow like ~ t^, as expected from equation 13.561 



3.2.4 Illustration of the modes of libration in a two planet 
system 

We now go on to consider two planet systems. As an illustrative example we discuss 
the GJ876 system (see table 13. We begin by considering the evolution of the system 
without stochastic forcing in order to characterise the modes of libration of the resonant 
angles and other orbital parameters. In particular, we identify the fast and slow modes 
discussed in section [3.1.21 

The time evolution of the resonant angles and the eccentricities of the unperturbed 
system is plotted in figure 13. 4[ Clearly visible are the slow and fast oscillation modes. 
The fast mode, which is seen to have a period of about 1.4 years, dominates the libra- 
tions of 62 and 0i while also being present in those of 02- On the other hand the slow 
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Figure 3.2: Growth of the semi-major axis Aa, the periastron Atu and the eccentricity 
Ae for a single planet as a function of time. The initial orbital parameters of the planet 
are taken to be those of GJ876 b (see table [XT]) . One hundred simulations starting 
from the same initial conditions are averaged for each value of the diffusion coefficient. 
The diffusion coefficient D is 8.2 ■ 10~^cm^/s^ for the upper and 8.2 ■ 10~^cm^/s^ for 
the lower curves, respectively. Solid lines correspond to the analytic predictions, i.e., 
these are not fits. Dashed lines correspond to the analytic predictions for uncorrelated, 
random kicks. 
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Figure 3.3: Growth of the difference in longitude AA as a function of time for a single 
planet. The initial orbital parameters of the planet are taken to be those of GJ876 b 
(see table 13. ip . One hundred simulations starting from the same initial conditions 
are averaged. Here, the diffusion coefficient is D = 8.2 ■ 10~^cm^/s^. The solid line 
corresponds to the analytic prediction, given by equation 13.561 



mode, which is seen to have a period of about 10 years, dominates the librations of C, 
while also being present in those of e\ and 02- 

We emphasise the fact that the eccentricities of the two planets participate in the 
librations, and are therefore not constant. In particular, the eccentricity of GJ876 b 
oscillates around a mean value of 0.03 with an amplitude Ae ~ 0.01 — 0.02. This 
behaviour involving the attainment of smaller values of the eccentricity has important 
consequences for stochastic evolution as discussed above (section I3.1.3P and also below. 

We remark that similar behaviour occurs for all the systems we have studied, these 
having a wide range of eccentricities and planet masses. The fast and slow mode periods 
scale with the planet masses as given by equations 13.311 and 13.331 respectively. These 
indicate, as confirmed by our simulations, that if both masses are reduced by a factor 
A then the period of the fast mode scales as and the period of the slow mode scales 
as A. 

In figure 13.51 we plot the librations in h and k variables. The left plots show one 
period of the slow mode libration, corresponding to one complete circle around the 
origin. Super-imposed on that are slow mode librations with higher frequency. The 
frequency ratio is ~ 6, thus the hexagon shape. The even higher frequency modes 
correspond to the orbital period itself. The right plots show many libration periods. 



3.3 Two planets with stochastic forcing 

We now consider systems of two planets with stochastic forcing. For simplicity we 
begin by applying additional forces only to the outer planet. We have found that 
adding the same form of forcing to the inner planet tends to speed up the evolution 
by approximately a factor of two without changing qualitative details. For illustrative 
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Figure 3.4: Time evolution of the resonant angles and the eccentricities in the system 
GJ876 without turbulent forcing. The dominance of the fast mode with period ~ 
1.4 years in the oscillations of 0i, and the dominance of the slow mode with period 
~ 10 years in ( can be clearly seen. 



purposes we again start with the GJ876 system and consider two diffusion coefficients 
D = 8.2 ■ 10~^cm^/s^ and D = 8.2 ■ 10~^cm^/s^. In all of our simulations a constant 
value of D = 2{Fi)^Tc is used by maintaining both parameters (Fj) and Tc constant. 
For the cases considered here, there is little orbital migration so that this effect can be 
ignored. 

The time evolution of the eccentricities is plotted for two runs in the first and 
third plot in figure 13.61 We see fast oscillations superimposed on a random walk. The 
amplitude of the oscillations, as well as the mean value, change with time. 

Our simple analytic model assumes slowly changing background eccentricities and 
semi-major axes. Accordingly, it does not incorporate the oscillations of the eccentricity 
due to the resonant interaction of the planets. In order to make a comparison, we 
perform a time average over many periods to get smoothed quantities whose behaviour 
we can compare with that expected from equations 13.411 - [3. 45[ When this procedure is 
followed (see second and forth plot in figure [32]), the evolution is in reasonable accord 
with that expected from the analytic model provided that allowance is made for the 
importance of small values of ei in determining the growth of the libration amplitude 
of the angle between the apsidal lines of the two planets (see equation I3.45p . The 
presence of this feature results in the behaviour of the libration amplitude being more 
complex than that implied by a process governed by a simple diffusive random walk. 
This is discussed in the following sections. 
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Figure 3.5: Librations in the GJ876 system, visualised using h and k variables. The 
initial conditions are the same as in figure 13.41 The left plots show the evolution during 
one slow mode period, which is equivalent to about 6 fast mode periods. The right 
plots show the evolution over many libration periods. 

3.3.1 Disruption of a resonance in stages 

Systems with mean motion commensurabilities can be in many different configurations. 
Here we describe the important evolutionary stages as they appear in a stochastically 
forced system starting from a situation in which all the resonant angles show small 
amplitude libration. Observations of GJ876 suggest that the system is currently in 
such a state with the ratio of the orbital periods P1/P2 oscillating about a mean value 
of 2. 

Attainment of circulation of the angle between the apsidal lines 

When the initial eccentricity of the outer planet, ei, is small, the excitation of the 
libration amplitudes of the resonant angles readily brings about a situation where ei 
attains even smaller values. This then causes the periastron difference C, to undergo 
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Figure 3.6: Time evolution of the eccentricities in the GJ876 system with turbulent 
forcing included. The diffusion coefficient is D = 8.2 ■ lO^^cm^/s^ for the two upper 
plots and D — S.2 • 10~^cm^/s^ for the two lower plots. The first (uppermost) and the 
third plots each show a single run. The second and fourth (lowermost) plots show the 
time averaged eccentricities for six different runs. The averaging interval is 1000 years, 
corresponding to ~ 100 slow mode periods. 



51 



1 




time [years] 

Figure 3.7: Growth of the difference in the resonance angle 0i as a function of time. 
The initial orbital parameters are those of GJ876 (see table ISTT]) . Here, the diffusion 
coefficient is -D = 8.2- 10~^cm^/s^. The solid line corresponds to the analytic prediction, 
given by equation 13.661 

large oscillations and eventually reach circulation (see equation 13. 45p . Should stochastic 
forces cause ei to reach zero, on account of the coordinate singularity, zui becomes 
undefined. Note that this occurs without a large physical perturbation to the system. 
Thus, the occurrence of this phenomenon does not imply the system ceases to be in a 
commensurability. 

In this context we note that the eccentricity of GJ876 b initially is such that ei ~ 0.03 
with values ei ~ 0.01 often being attained during libration cycles. Thus only a small 
change may cause the above situation to occur. In all cases that we have considered, 
we find that ( enters circulation prior to the fast angle 0i, which may remain librating 
until that too is driven into circulation. 

Attainment of circulation of the fast angle 

Both before and after ( enters circulation, stochastic forces increase the libration am- 
plitude of the fast mode 0i. This mode dominates both the librations of the resonant 
angle (pi and the semi-major axes. Eventually 0i starts circulation. Shortly afterwards, 
commensurability is lost and P1/P2 starts to undergo a random walk with a centre that 
drifts away from 2. Note that it is possible for some realisations to re-enter commen- 
surability. For systems with the masses of the observed GJ876 system, the most likely 
outcome is a scattering event that causes complete disruption of the system. 

Because the eccentricities do not play an important role in the attainment of cir- 
culation of the fast angle, we can use equation 13.661 to calculate the growth rate. The 
root mean square values of A0i as a function of time, taken from both, analytical and 
numerical results, are plotted in figure 1X71 Here, we use D = 8.2 ■ 10~^cm^ /s^. Note 
that as mentioned above, ( is circulating first (in the presented case this happens after 
~ 1000 years). However, this does not influence the growth of A(f)i. 
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3.3.2 A numerical illustration 



In order to illustrate the evolutionary sequence described above we plot results for 
two realisations of the evolution of the GJ876 system in figure 13.81 For these runs 
we adopted the diffusion coefficient D = 0.42cm^/s^. Note that increasing D only 
decreases the evolutionary time which has been found, both analytically and numerically 
to be oc 1/D (see below). 

The times at which the transition from libration to circulation occurs for both the 
slow and fast angles, are indicated by the vertical lines in figure 13.81 

Several of the features discussed above and in section 13.1.31 can be seen in figure 
13.81 In particular the tendency for the occurrence of very small values of ei to be 
associated with transitions to circulation of ( can be seen for the realisation plotted 
in the lower panels at around t = 80 years and t = 500 years. Such episodes always 
seem to occur when the libration amplitude of the fast angle 0i is relatively boosted, 
indicating that this plays a role in boosting the slow angle. If the duration for which 
ei attains small values is small and (pi recovers small amplitude librations, the slow 
angle returns to circulation. Thus, the attainment of long period circulation for the 
slow angle is related to the diffusion time for 01. We also see from figure ESI that the 
angle 02? which has a large contribution from the slow mode, behaves in the same way 
as ( as far as libration/circularisation is concerned. We have verified by considering the 
results from the simulations of GJ876 LM HE, which started with a larger value of ei, 
that, as expected, the attainment of circulation of the slow angle takes relatively longer 
in this case, the time approaching more closely the time when 0i attains circulation. 
Also as expected, the time when 0i attains circulation is not affected by the change in 
ei. 

3.3.3 Dependence on the diffusion coefficient 

We now consider the stability of the systems listed in table 13.11 as a function of D. 
These systems have a variety of masses and orbital eccentricities. In particular, in view 
of the complex interaction between the resonant angles discussed above, we wish to 
investigate whether the mean amplitude growth at a given time is indeed proportional 
to D. As also mentioned above, the value of the diffusion parameter D that should be 
adopted, is very uncertain. We have therefore considered values of D ranging over five 
orders of magnitude. The correlation time Tc is always taken to be given by Tc = 0.5Q~^ 
while the RMS value of the force is changed. In order to determine the lifetime of a 
resonant angle, we monitor whether it is librating or circulating. Numerically, libration 
is defined to cease when the angle is first seen to reach absolute values larger than 2. 
We note that the angle can in general regain small values afterwards. However, this is 
a transient effect and changes the lifetime by no more than a factor of ~ 2 in all our 
simulations. 

In this context we consider the fast angle (pi and the slow angle (. As we saw above, 
the latter can be replaced by (p2 which exhibits the same behaviour. Because (pi is the 
last to start circulating, the resonance is defined to be broken at that point. 

Equations 13.661 and 13.651 estimate the spreading of the resonant angles as a function 
of time. We plot both the numerical and analytical results in figure 13.91 To remove 
statistical fluctuations and obtain a mean spreading time, the numerical values for a 
particular value of D were obtained by averaging over 60 realisations. 
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Figure 3.8: Time evolution of the resonant angles, the period ratio P2/P1 and the 
eccentricity, ei, in the GJ876 system with stochastic forcing corresponding to D = 
0.42cm^/s'^. The vertical lines indicate when the angles enter circulation for a prolonged 
period. The realisation illustrated in the lower panel scatters shortly after 0i goes into 
circularisation. 
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Figure 3.9: Average time until circulation of the resonant angles 0i (top) and 02 (bot- 
tom) in the systems GJ876 (indicated with +) , GJ876LM (indicated with x), GJ876SE 
(indicated with ★) and G J876E (indicated with □) as a function of the stochastic forcing 
diffusion coefficient D. In the case of the upper panel, the analytic curves explained 
in the text, are from top to bottom for the GJ876, GJ876LM, GJ876SE and GJ876E 
systems, respectively. In the lower panel the analytic curves, as explained in the text, 
are from top to bottom for the GJ876 (solid curve and top dashed curve), GJ876LM, 
GJ876SE and GJ876E systems, respectively. 
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From figure 13.91 it is apparent that the evolutionary times scale is oc 1/Z^ for 
varying by many orders of magnitude. The analytic estimates for the libration survival 
times of 0i, dominated by the fast mode, obtained using equation 13.661 adopting the 
initial orbital elements and with the fast libration frequency determined from the simu- 
lations, are plotted in the upper panel of figure 13.91 These are in good agreement with 
the numerical results. However, the estimate for 02 based on equation 13.651 using the 
initial value of ei overestimates the lifetime by a factor of at least ~ 40 (see solid line in 
the lower plot of figure l3^ . Furthermore equation 13.651 has no dependence on planet 
mass which is in clear conflict with the numerical results. As discussed above, this is 
due to the temporary attainment of small eccentricities. This causes the disruption 
of the libration of ( earlier than predicted assuming ei being constant. In fact this 
disruption occurs at times that can be up to 10 times shorter than those required to 
disrupt 01. We estimate the lifetime of librations of 02 by calculating the time 0i needs 
to reach values close to 1 (see dashed lines in the lower plot of figure [3l9|) . As explained 
above in section I3.1.3[ large amplitude variations of 0i are expected to couple to and 
excite the slow mode. Variations induced in the eccentricity ei allow ei to reach zero 
and we lose libration of 02. These simple estimates are in good agreement with the 
numerical results and accordingly support the idea that 0i and 02 are indeed coupled 
in the non-linear regime. For low mass planets this simple analytic prediction underes- 
timates the lifetime of 02 by a factor of ~ 2, suggesting that the coupling between the 
two modes depends shghtly on the planet masses. 

To confirm our understanding of these processes, we have repeated the calculations 
with systems that are in the same state as the system discussed above but have higher 
eccentricities (see GJ876 LM HE and GJ876 SE HE in table El]). The lifetime of 
the librating state of 0i is unchanged, as this does not depend significantly on the 
eccentricity. However, the lifetime of the librating state of 02 is a factor of 2 — 3 longer. 
This is due to the fact that a larger excitation of 02 is needed to make ei reach small 
values in these simulations. 

The random walk description breaks down completely when the anticipated disrup- 
tion time becomes shorter than the libration period. This is expected to happen for 
small enough masses as can be verified from figure 13.91 It is because the disruption 
time decreases linearly with the planet masses while the libration period increases as 
the square root of the planet masses. Then we cannot average over many libration 
periods. This situation is apparent in figure 13.91 for very short disruption times of the 
order 100-1000 years where survival times cease to vary linearly with 1/D. 



3.4 Formation of HD128311 



Having understood the effects of stochastic forces on resonances, we now go on to 
discuss the application of these ideas to the formation of the orbital configuration of 
the planetary system HD128311. This is (with 99% confidence) in a 2:1 mean motion 
resonance, with the angle 0i libratin g and the angle ( circulating so that there is no 
apsidal co-rotation f lVogt et al.ll2005l ). Ho wever, the o rbital parameters are not well 
constrained. The original Keplerian fit by IVogt et al.l (|2005[ ) is such that the planets 
undergo a close encounter after only 2000 years. The values in table 13.11 have been 
obtained from a fit to the observational data that includes interactions between the 
planets. The values quoted have large error bars. For example the eccentricities ei and 
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62 have an uncertainty of 33% and 21%, respectively. Although the best fit doesn't 
manifest apsidal co-rotation, the system could undergo large amphtude librations and 
still be stable. 

According to iLee fc Peald (120021 ) the planets should have apsidal co-rotation if the 
commensurability was formed by the two planets undergoing sufficiently slow convergent 
inward migration, and if they then both migrated inward s significantly while maintain- 
i ng th e commensurability (see also ISnellgrove et al.ll200ll ). Accordingly ISandor &: Kiev 
(120061 ) suggested a possible formation scenario with inward migration as described 
above, but with an additional perturbing event, such as a close encounter with an 
additional third planet occurring after the halting of the inward migration. This per- 
turbation is needed to alter the behaviour of (, so that it undergoes circulation rather 
than libration and thus producing orbital parameters similar to the observed ones. 

We showed above that stochastic forcing possibly resulting from turbulence driven by 
the MRI readily produces systems with commensurabilities without apsidal co-rotation 
if the eccentricities are not too large. This suggests that a scenario that forms the com- 
mensurability through disc-induced inward convergent migration might readily produce 
commensurable systems without apsidal co-rotation if stochastic forcing is included. 
Such scenarios are investigated in this section. 

We present simulations of the formation of HD128311 that include migration and 
stochastic forcing due to turbulence but do not invoke artificial perturbation events 
involving additional planets. We find that model systems with orbital parameters re- 
sembling th e observed ones are re adily produced and are better matches than those 
provided by lSandor fc Kleyl ( 120061 ). The planets in this system are of the order of sev- 
eral Jupiter masses and the eccentricity of one planet can get very small during one 
libration period. The observed orbital parameters are given in table 13.11 and their time 
evolution (due to planet-planet interactions) is plotted on the right hand side of figures 
13101 KW and Km The mass of the star is 0.8 Mq. 



3.4.1 Migration forces 

We incorporate the non-conservative forces exerted by a proto-stellar disc that lead to 
inward migration by following the approach outlined in appendix [Bl The migration 
process is characterised by migration and eccentricity damping times cales, t„, = a/ a 



and Te = e/e, respecti vely. This procedure is now widely used (e.g. ISandor fc Kley 



20061 : ICrida et al.ll2008[ ). We also use the ratio of the above timescales, K = \Ta \ /re, as 
defined in chapter [2J This ratio determines the eccentricities in the state of self-similar 
inward migration of the commensurable system that is attained after large times. 

In this work we allow the planets to form a commensurability through convergent 
inward migration with stochastic forcing included. The disc is then removed through 
having its surface density reduced to zero on a 2000 year timescale. This reduces both 
the migration forces and the stochastic forces simultaneously. This procedure is very 
similar to that adopted by the above authors but we have included stochastic forcing 
and removed the disc on longer, more realistic, time-scales. 

The stochastic forces need to have the right balance with respect to the migration 
rate. We have found that inward migration imparts stability to the resonant system. 
If the migration rate is too fast relative to the stochastic forcing the migration keeps 
down the libration amplitudes and we do not get circulation. On the other hand 



57 




5000 10000 15000 20000 ' C 4000 

t (years) t (years) 



Figure 3.10: The plots on the left hand side show a possible formation scenario for the 
system HD 128311 including turbulence and migration. We plot the observed system on 
the right hand side as a comparison (see table [XT] and text). The plots show, from top 
to bottom, the resonant angles (j)i,(f)2,C^ the eccentricities 61,62 and the period ratio 
Pi/P2- Resonance capturing occurs after 2000 years. 



large eccentricity damping favours broken apsidal co-rotation. This might look counter 
intuitive at first sight, but remember that the diffusion of ( depends on 1/6^ and 1/63. 

Due to the stochastic nature of the problem, it is hard to present a continuum of 
solutions so we restrict ourselves to the discussion of three representative examples. 
However, we comment that we are able to obtain similar end states for a wide range of 
migration parameters. 

As above, the parameters (Fj) and Tc are kept constant, so maintaining D constant, 
for the duration of the simulation, with Tc being determined for the initial location of 
the outer planet. As discussed above, these values may scale with the radius of the 
planet and we thus expect them to change during migration. However, the semi-major 
axis of the outer planet changes only by ~ 30 % during the simulation. Consequently, 
we ignore this effect. 

3.4.2 Model 1 

The planets are initially on circular orbits at radii ai = 4 AU and 02 = 2 AU. Stochastic 
forcing is applied to the outer planet only, with a diffusion coefficient of D = 6.4 • 
10~^cm^/s^. Note that although this value is 64 times larger than the scaling given by 
equation I3.71[ corresponding to a force that is 8 times larger than the simple estimate 
given in section l3.2.H corresponding to smaller reduction factors resulting from a gap or 
dead zone, it is too small for the angle 0i to be brought to circulation during our runs. 
The results given above and summarised in figure 13. 9[ and further tests, indicate that 
similar results would be obtained if D is reduced, but on a longer time-scale oc 1/D, 
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provided the migration rate is also ultimately reduced appropriately so that the system 
can survive for long enough to enable ( to be driven into circulation. 

The outer planet is made to migrate inwards on a timescale r^^ = 8000 years. The 
inner planet migrates slowly outwards on a timescale of Tag = —20000 years. This results 
in convergent migration. For both planets we use an eccentricity to semi-major axis 
damping ratio of K = 8. The resu l ting e volutionary timescales are significantly larger 
than those used by ISandor Sz KlevI fcooeh and more easily justified by hydrodynamical 
simulations. 

The time evolution is shown in the left plot of figure 13.101 After resonance capturing 
all angles are either initially librating or on the border between libration and circulation. 
The slow mode has a period of ~ 700 years, whereas the fast mode period is 20 times 
shorter. 

It can be seen from figure [3?T0] and other similar figures below, that while migration 
continues, libration amplitudes tend to be controlled apart from when ei becomes either 
zero or very close to zero. Then, due to stochastic forcing, 02 and ( start circulating. 
Subsequently, libration is recovered over time intervals for which ei does not attain very 
small values but eventually additional stochastic forcing together with the repeated 
attainment of small values for ei causes 02 to remain circulating for the remainder of 
the simulation. After 13000 years, both the forces due to migration and turbulence are 
reduced smoothly on a timescale of 2000 years. The result is a stable configuration that 
resembles the observed system very well. 



3.4.3 Model 2 

For this simulation, we lengthen the migration timescales by a factor of 2 to show that 
the results are generic. We also decrease the value of K to 5.5. All other parameters 
are the same as for model 1. This results in larger eccentricities and the final state 
better resembled our representation of the observed system. However, it should be kept 
in mind that the eccentricities are not well constrained by the observations. 

The time evolution is plotted in figure I3TTT1 The resonant angles librate immediately 
after the capture into resonance, as predicted for sufficiently slow migration. However, 
in the same way as for model 1 described above, stochastic forces make 02 circulate soon 
afterwards. Once the migration forces and stochastic forces are removed between 20000 
and 22000 years, the system stays in a stable configuration with no apsidal co-rotation. 



3.4.4 Model 3 

All parameters are the same as for model 2. Accordingly, model 3 represents another 
statistical realisation of that case. The time evolution is plotted in figure 1X121 The final 
state is very close to the boundary between libration and circulation of C- As discussed 
above, large uncertainties in the orbital parameters mean that the actual system could 
be in such a state. 

Due to stochastic forces, which may result from local turbulence, we are able to 
generate a broad spectrum of model systems. Some of them undergo a strong scattering 
during the migration process. However, the surviving systems resemble the observed 
system very well. Although the orbital parameters are not well constrained yet, we 
showed that values in the right range are naturally produced by a turbulent disc. 
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Figure 3.11: The plots on the left hand side show another possible formation scenario 
for the system HD128311. Again, we plot the observed system on the right hand side 
as a comparison (see table 13.11 and text). The migration timescales in this run are 
Ta,i = 16000, Ta,2 = —40000 and K = 5.5. Note that the eccentricities are larger 
compared to figure |3TT0] because K is smaller. 
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Figure 3.12: These plots show another possible formation scenario of HD128311. The 
damping parameters are the same as in figure 13.111 
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3.5 Conclusions 



In this chapter, we have presented a self-consistent analytical model of both a single 
planet, and two planets in a mean motion resonance, being subject to external stochastic 
forcing. The stochastic forces could result from MRI-driven turbulence within the proto- 
planetary disc but our treatment is equally applicable to any other source, for example 
from a gravitoturbulent disc such as Saturn's A-ring (see chapter H]). 

We considered the evolution of a stochastically forced two planetary system that is 
initially deep inside a MMR (i.e., the two independent resonant angles librate with small 
amplitude). Stochastic forces cause libration amplitudes to increase in the mean with 
time until all resonant angles are eventually driven into circulation, at which point the 
commensurability is lost. Often a strong scattering occurs soon afterwards for systems 
composed of planets in the Jovian mass range. 

We isolated a fast libration mode which is associated with oscillations of the semi- 
major axes and a slow libration mode which is mostly associated with the motion of 
the angle between the apsidal lines of the two planets. These modes respond differently 
to stochastic forcing, the slow mode being more readily converted to circulation than 
the fast mode. This slow mode is sensitive to the attainment of small eccentricities, 
which cause rapid transitions between libration and circulation. The amplitude of the 
fast mode grows more regularly in the mean, with the square of the libration amplitude 
in most cases increasing linearly with time and being proportional to the diffusion 
coefficient D. Of course this discussion is simplified and there are limitations. For 
example if the total mass of the system is reduced, the disruption time eventually 
becomes comparable to the libration period. In that case the averaging process that we 
use in the derivation is no longer valid and the lifetime no longer scales as 1/D. 

The analytic model was compared to numerical simulations which incorporate stochas- 
tic forces. Those forces, parametrised by the mean square values of each component 
in cylindrical coordinates and the auto-correlation time, are applied in a continuous 
manner giving results that could be directly compared with the analytic model. The 
simulations are in broad agreement with analytic predictions and we presented illustra- 
tive examples of the disruption process. We performed the simulations for a large range 
of diffusion parameters, planet masses and initial eccentricities to verify the scaling law 
for the commensurability disruption time summarised hereafter. 

To summarise our results, recalling that the slow angle is driven into circulation 
before the fast angle, so that the ultimate lifetime is determined by the time taken for 
the fast angle to achieve circulation, we determine the lifetime, t/, using equation I3.66[ 
setting t = tf, (A0i)^ = (A</)i)g = 4, together with 7^ = p = 1, we obtain 

= ^ '^^ 3.73 
^ 36^ ^ ' 

We showed above that this gives good agreement with our numerical results. Using 
D = 2{F^)tc, we can express this result in terms of the relative magnitude of the 
stochastic forcing in the form 

2^4\ //AJ,^2\ /QK,.,.. /7rrz\'^ 



tj = 2.4 X 10-^ ( ^Sv ) ( ) ■ ( ) —^1, (3.74) 



/(A0i)g\ f8.5uif,/q^\ q 



where Pi is the orbital period of the outer planet. Here the first quantity in brackets 
represents the ratio of the square of the central force per unit mass to the mean square 
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stochastic force per unit mass acting on the outer planet. The other quantities in 
brackets, scaled to the GJ876 system are expected to be unity, while the last factor 
q/qcj is the ratio of the total mass ratio of the system to the same quantity for GJ876. 
Here, it is assumed that the two planets in the system have comparable masses. 

From equation 13. 74[ we see that a non-migrating system such as GJ876 could survive 
in resonance for t/ ~ 10^ years if the stochastic force amplitude is ~ 10^^ times the 
central force. This expression enables scaling to other systems at other disc locations 
for other stochastic forcing amplitudes. Inference of survival probabilities for particular 
systems depends on many uncertain aspects, such as the proto-planetary disc model 
and the strength of the turbulence. However, the mass ratio dependence in equation 
13.741 indicates that survival is favoured for more massive systems. At the present time 
the number of observed resonant systems is too small for definitive conclusions to be 
made. However, the fact that several systems exhibiting commensur abilities have been 
observed indicates that resonances are not always completely disrupted by stochastic 
forces due to turbulence, but rather may be modified as in our study of HD128311. 

The most likely orbital configuration of the HD128311 system is such that the fast 
mode librates with the slow mode being near the borderline between libration and 
circulation. We found that such a configuration was readily produced in a scenario 
in which the commensurability was formed through a temporary period of convergent 
migration, with the addition of stochastic forcing. During a migration phase moderate 
adiabatic invariance applied to the libration modes together with eccentricity damping 
leads to increased stabilisation and a longer lifetime for the resonance. However, the 
time evolution of the eccentricity, and in particular the attainment of small values, plays 
an important role in causing the slow mode to circulate. This corresponds to the loss 
of apsidal co-rotation. Thus we expect that a large eccentricity damping rate does not 
necessarily stabilise the apsidal co-rotation of a system. Additional simulations have 
shown that this is lost more easily for large damping rates {K ^ 10). On the other 
hand, it should be noted that this has less significance when one of the orbits becomes 
nearly circular. 

Further observations of extra-solar planetary systems leading to better statistics 
may lead to an improved situation for assessing the role of stochastic forcing. 
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Chapter 



Stochastic migration of small bodies in 

Saturn's rings 



Have not the small particles of bodies certain powers, 
virtues, or forces, by which they act at a distance, not 
only upon the rays of light for reflecting, refracting, 
and inflecting them, but also upon one another for 
producing a great part of the Phenomena of nature? 

Isaac Newton, Optiks, Query 31, 1729 



In the previous chapter, we discussed the stochastic migration of planets, embedded 
in a circum-stellar disc. A very similar process also occurs in the circum-planetary 
disc known as Saturn's rings. However, both the length and time scales involved 
are three orders of magnitude smaller. The bodies of interest are small moonlets. 
These 20 m — 100 m sized bodies have been predic t ed both analytic ally and numerically 



(ISpahn fc Sremcevid l2000l : ISremcevic et al.l l2002l : ISeifi et al.l l2005l ) . They create pro- 



peller shaped structures due to their disturbance of the rings a nd have been observed - 



only recently by the Cassini spacecraft in the A and B ring (ITiscareno et al.l 12006 



20081 ). 



They can migrate within the rings, similar to proto-planets which migrate in a 
proto-stellar disc. Depending on the disc properties and the moonlet size, this can 
happen in either a smooth or in a stochastic (random walk) fashion. We refer to those 
migration regimes as type I and type IV, re spectively, iri analogy to the terminology in 
disc-planet Weractions (see section 03. Icricia ct al I H) 

showed that there is a 

laminar type I regime that might be important on very long timescales. This migration 
is qualitatively different to the migration in a pressure supported gas disc. However, 
the migration of moonlets in the A ring is generally dominated by type IV migration, 
at least on short timescales. 

We here study their dynamical evolution in the type IV regime analytically, including 
important effects such as collisions which were not present in the discussion in chapter 
Ei We also perform realistic three dimensional collisional N-body simulations with up to 
a quarter of a million particles. A new set of pseudo shear periodic boundary conditions 
is used which reduces the computational costs by an order of magnitude compared to 
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previous studies. The numerical simulations confirm our analytic estimates to within a 
factor of two. 

For low ring surface densities, the main effects on the evolution of the eccentric- 
ity and the semi-major axis are found to be due to collisions and the gravitational 
interaction with particles in the vicinity of the moonlet. For large surface densities, 
the gravitational interaction with self-gravitating wakes becomes important. These 
stochastic forces are very similar to what has been discussed earlier, in chapter |3l 

On short timescales the evolution is always dominated by stochastic effects caused 
by collisions and gravitational interaction with self-gravitating ring particles. These 
result in a random walk of the moonlet's semi-major axis. The eccentricity of the 
moonlet quickly reaches an equilibrium value due to collisional damping. The average 
change in semi-major axis of the moonlet after 100 orbital periods is 10-lOOm. This 
translates to an offset in the azimuthal direction of several hundred kilometres. We 
expect that such a shift is easily observable. 

We start be reviewing the basic equations governing the moonlet and the ring par- 
ticles in a shearing box approximation in section 14.11 Then, in section \4.2\ we estimate 
the eccentricity damping timescale due to ring particles colliding with the moonlet and 
ring particles on horseshoe orbits as well as the effect of particles on circulating orbits. 
In section we estimate the excitation of the moonlet eccentricity caused by stochas- 
tic particle collisions and gravitational interactions with ring particles. This enables us 
to derive an analytic estimate of the equilibrium eccentricity. 

In section 1131 we discuss and evaluate processes, such as collisions and gravitational 
interactions with ring particles and self-gravitating clumps, that lead to a random walk 
in the semi-major axis of the moonlet. 

We describe our numerical setup and the initial conditions used in section 14.51 The 
results are presented in section H^l All analytic estimates are confirmed both in terms of 
qualitative trends and quantitatively to within a factor of about two in all simulations 
that we performed. We also discuss the long term evolution of the longitude of the 
moonlet and its observability, before summarising our results in section 14.71 

4.1 Basic equations governing the moonlet and 
ring particles 

4.1.1 Epicyclic motion 

We adopt a local right handed Cartesian coordinate system with its origin being in 
circular Keplerian orbit with semi-major axis a and rotating uniformly with angular 
velocity fl. This orbit coincides with that of the moonlet when it is assumed to be 
unperturbed by ring particles. The x axis coincides with the line joining the central 
object of mass Mp and the origin. The unit vector in the x direction, e^., points away 
from the central object. The unit vector in the y direction points in the direction of 
rotation and the unit vector in the z direction, Sz points in the vertical direction being 
normal to the disc mid-plane (see also appendix |D|) . 

In general, we shall consider a ring particle of mass mi interacting with the moonlet 
which has a much larger mass m2. A sketch of three possible types of particle trajectory 
in the vicinity of the moonlet is shown in figure 14.11 These correspond to three distinct 
regimes, a) denoting circulating orbits, b) and c) denoting orbits that result in a collision 
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Figure 4.1: Trajectories of ring particles in a frame centred on the moonlet. Particles 
accumulate near the moonlet and fill its Roche lobe. Particles on trajectories labelled a) 
are on circulating orbits. Particles on trajectories labelled b) collide with other particles 
in the moonlet's vicinity. Particles on trajectories labelled (c) collide with the moonlet 
directly. Particles on trajectories labelled (d) are on horseshoe orbits. Particles on 
trajectories labelled (e) leave the vicinity of the moonlet. 



with particles in the vicinity of the moonlet and directly with the moonlet respectively, 
and d) denoting horseshoe orbits. All these types of trajectory occur for particles that 
are initially in circular orbits, both interior and exterior to the moonlet, when at large 
distances from it. 

Approximating the gravitational force due to the central object by its first order 
Taylor expansion about the origin leads to Hill's equations (see equations ID.SIIDTToI) . 
governing the motion of a particle of mass mi of the form 

r = -2fie^ X r + 3fi^(r ■ e^)e^ - V^f/mi, (4.1) 

where r = (x, z) is position of a particle with mass mi and \l/ is the gravitational 
potential acting on the particle due to other objects of interest such as the moonlet. 
The square amplitude of the epicyclic motion S"^ can be defined through 

£^ = n-^x^ + {2n'^y + 3xy. (4.2) 

Note that neither the eccentricity e, nor the semi-major axis a are formally defined 
in the local coordinate system. However, in the absence of interaction with other 
masses (\E' = 0), £^ is conserved, and up to first order in the eccentricity we may make 
the identification S = ea. We recall the classical definition of the eccentricity e in a 
coordinate system centred on the central object 



w X s X w s 



n^a^ Isl 



(4.3) 



Here, s is the position vector (a, 0, 0) and w is the velocity vector of the particle 
relative to the mean shear associated with circular Keplerian motion as viewed in the 
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local coordinate system plus the circular Keplerian velocity corresponding to the orbital 
frequency Q of the origin. Thus 

w = ^x, y + ^rix + an, . (4.4) 

All eccentricities considered here are very small (~ 10~* - 10~'^) so that the differ- 
ence between the quantities defined through use of equations 14.21 and 14.31 is negligible, 
allowing us to adopt £^ as a measure of the eccentricity throughout this chapter. 

Let us define another quantity, A, that is also conserved for non interacting particle 
motion, \E' = 0, which is the x coordinate of the centre of the epicyclic motion and is 
given by 

A = 2Q-^y + 4:X. (4.5) 

We identify a change in ^ as a change in the semi-major axis a of the particle, again, 
under the assumption that the eccentricity is small. 



4.1.2 Two interacting particles 

We now consider the motion of two gravitationally interacting particles with position 
vectors ri = {xi,yi,zi) and r2 = (2:2, 1/2, 2^2) • Their corresponding masses are mi and 
7712, respectively. The governing equations of motion are 

ri = -2(^e^ X ri + 3fi^(ri-e^)e^- Vri^i2/mi (4.6) 
^2 = -2^6, X 1-2 + 3fi2(r2 ■ e^.)ex- + Vri^i2/m2, (4.7) 

where the interaction gravitational potential is \E'i2 = —Gmim2/ \ri — r2|. The position 
vector of the centre of mass of the two particles is given by 

77Jiri + 7772r2 

r = . (4.8) 

mi + 7/72 

The vector r also obeys equation 14.11 with = 0, which applies to an isolated particle. 
This is because the interaction potential does not affect the motion of the centre of 
mass. We also find it useful to define the vector 

£, = {Q-^Xi, 2Q-^yi + 3xi) i = 1,2. (4.9) 

Then, consistently with our earlier definition of we have Si = \£i\. The amplitude of 
the epicyclic motion of the centre of mass S is given by 

(7r7i + 7^72)^^^ = mlS'^ + m\£l + 2mim2£i£2 cos(0i2). (4-10) 

Here (/)i2 is the angle between £1 and £2- It is important to note that £ is conserved 
even if the two particles approach each other and become bound. This is as long as 
frictional forces are internal to the two particle system and do not affect the centre of 
mass motion. 
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4.2 Effects leading to damping of the eccentricity 
of the moonlet 



We begin by estimating the moonlet eccentricity damping rate associated with ring 
particles that either collide directly with the moonlet, particles in its vicinity, or interact 
only gravitationally with the moonlet. 

4.2.1 The contribution due to particles impacting the 
moonlet 

Particles impacting the moonlet in an eccentric orbit exchange momentum with it. Let 
us assume that a ring particle, identified with mi has zero epicyclic amplitude, so that 
Si = far away from the moonlet. The moonlet is identified with m2 and has an initial 
epicyclic amplitude 82- The epicyclic amplitude of the centre of mass is therefore 



where we have assumed that mi ^ 1712. 

The moonlet is assumed to be in a steady state in which there is no net accretion 
of ring particles. Therefore, for every particle that either collides directly with the 
moonlet or nearby particles bound to it (and so itself becomes temporarily bound to it), 
one particle must also escape from the moonlet. This happens primarily through slow 
leakage from locations close to the Li and L2 points such that most particles escape from 
the moonlet with almost zero velocity (as viewed from the centre of mass frame) and 
so do not change its orbital eccentricity. However, after an impacting particle becomes 
bound to the moonlet, conservation of the epicyclic amplitude associated with the 
centre of mass motion together with equation 14.111 imply that each impacting particle 
will reduce the eccentricity by a factor 1 — mi/m2. 

It is now an easy task to estimate the eccentricity damping timescale by determining 
the number of particle collisions per time unit with the moonlet or particles bound to 
it. To do that, a smooth window function Wh+c{b) is used, being unity for impact 
parameters b that always result in an impact with the moonlet or particles nearby that 
are bound to it, being zero for impact parameters that never result in an impact. 

To ensure that the impact band is well defined, we assume that the epicyclic am- 
plitude is smaller than the moonlet 's Hill radius, which is the case in all simulations 
discussed below. However, if this were not the case, one could simply calculate the 
collisions per time units by other means, for example by using the velocity dispersion 
and the geometrical cross section of the moonlet. 

The number of particles impacting the moonlet per time unit, dN/dt, is obtained 
by integrating over the impact parameter with the result that 



We note that allowing b to be negative enables impacts from both sides of the moonlet 
to be taken into account. 




(4.11) 




(4.12) 
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Therefore, after using equation 14. Ill we find that the rate of change of the moonlet's 
eccentricity 62, or equivalently of it's amphtude of epicychc motion £2, is given by 

d£2 82 £2 r 3 



dt ''"e, collisions ITI2 



/oo o 
Wb+c{b)db, (4.13) 
-00 ^ 



where Te,coiiisions defines the circularisation time arising from colhsions with the moon- 
let. We remark that the natural unit for b is the Hill radius of the moonlet, = 
(m2/(3Mp))^/^a, so that the dimensional scaling for Te^coiusions is given by 

Cm.ions « GJ:r-^'Q^' (4.14) 

which we find to also apply to all the processes for modifying the moonlet's eccentricity 
discussed below. If we assume that PFt+cd^l) can be approximated by a box function, 
being unity in the interval [l.brfj, 2.5r/^] and zero elsewhere, we get 

Cms.ons = 2.0 G J: rj,' n-\ (4.15) 

4.2.2 Eccentricity damping due the interaction of the 
moonlet with particles on horseshoe orbits 

The eccentricity of the moonlet manifests itself in a small oscillation of the moonlet 
about the origin. Primarily ring particles on horseshoe orbits will respond to that 
oscillation and damp it. This is because only those particles on horseshoe orbits spend 
enough time in the vicinity of the moonlet, i.e. many epicyclic periods. 

In appendix [Cl we calculate the amplitude of epicyclic motion £if (or equivalently 
the eccentricity ei/) that is induced in a single ring particle in a horseshoe orbit un- 
dergoing a close approach to a moonlet which is assumed to be in an eccentric orbit. 
In order to calculate the circularisation time, we have to consider all relevant impact 
parameters. Note that each particle encounter with the moonlet is conservative and 
is such that for each particle, the Jacobi constant, applicable when the moonlet is in 
circular orbit, is increased by an amount by the action of the perturbing 

force, associated with the eccentricity of the moonlet, as the particle passes by. Because 
the Jacobi constant, or energy in the rotating frame, for the moonlet and the particle 
together is conserved, the square of the epicyclic amplitude associated with the moon- 
let alone changes by Sfjrmi/m2. Accordingly the change in the amplitude of epicyclic 
motion of the moonlet 82, consequent on inducing the amplitude of epicyclic motion 
Sif in the horseshoe particle, is given by 

Note that this is different compared to equation 14.111 Here, we are dealing with a 
second order effect. First, the eccentric moonlet excites eccentricity in a ring particle. 
Second, because the total epicyclic motion is conserved, the epicyclic motion of the 
moonlet is reduced. 

Integrating over the impact parameters associated with ring particles and taking into 
account particles streaming by the moonlet from both directions by allowing negative 
impact parameters gives 



d£i 



dt 



horseshoe 



3E£',fn\b\W,{b)db 2£l 

= , (4.17) 

■^^2 ^e, horseshoe 
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where Tg^horseshoe IS the circularisation time and, as above, we have inserted a window 
function, which is unity on impact parameters that lead to horseshoe orbits, otherwise 
being zero. Using £if given by equation IC.181 we obtain 



_i _ 9 /Sri^] 

'"e, horseshoe 



128 (^) rj'"^'"^' {i^VTHY") dv. (4.18) 

For a sharp cutoff of VF^dfol) at bm = l-Sr/^ we find the value of the integral in the 
above to be 2.84. Thus, in this case we get 

^iseshoe = 0.13 G S r^' n-\ (4.19) 

However, note that this value is sensitive to the value of bm adopted. For bm = l.25rH, 
Te is a factor of 4.25 larger. 



4.2.3 The effects of circulating particles 



The effect of the response of circulating particles to the gravitational pe rturbation of the 

moonlet on the moonlet's eccentrici t y can be estimated from the work of lGoldreich &: Tremaine 
( 1l980l ) and iGoldreich fc Tremaind (119821 ) . These authors considered a ring separated 
from a satelhte such that co-orbital effects were not considered. Thus, their expressions 
may be applied to estimate effects due to circulating particles. However, we exclude 
their corotation torques as they are determined by the ring edges and are absent in 
a local model with uniform azimuthally averaged surface density. Equivalently, one 
may simply assume that the corotation torques are saturated. When this is done only 
Lindblad torques act on the moonlet. These tend to excite the moonlet's eccentricity 
rather than damp it. 



We replace the ring mass in equation 70 of IGoldreich &: Tremaind ( 119821 ) by the 
integral over the impact parameter and switch to our notation to obtain 



dt 



9.55/ m2T.G'n-'\b\-'> £^Wa{h)db 



(4.20) 



where Wa{b) is the appropriate window function for circulating particles. Assuming a 
sharp cutoff of M4(|6|) at 6^, we can evaluate the integral in equation 14.201 to get 



l_ d£2 

82 dt 



1 



0.183 G S 



(4.21) 



where we have adopted bm = 2.5rj|/, consistently with the simulation results presented 
below. We see that, although Te^circ < corresponds to growth rather than damping 
of the eccentricity, it scales in the same manner as the circularisation times in sections 
14.2.11 and 14.2.21 scale (see equation I4.14p . 

This timescale and the timescale associated with particles on horseshoe orbits, 
7"e,horseshoe5 are significantly larger than the timescale associated with collisions, given in 
equation 14.151 We can therefore ignore those effects for most of the discussion in this 
chapter. 
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4.3 Processes leading to the excitation of the 
eccentricity of the moonlet 



In section 14.21 we assumed that the moonlet had a small eccentricity but neglected the 
initial eccentricity of the impacting ring particles. When this is included, collisions 
and gravitational interactions of ring particles with the moonlet may also excite its 
eccentricity. 



4.3.1 Collisional eccentricity excitation 

To see this, assume that the moonlet initially has zero eccentricity, or equivalently no 
epicyclic motion, but the ring particles do. 

As above we consider the conservation of the epicyclic motion of the centre of mass 
in order to connect the amplitude of the final epicyclic motion of the combined moonlet 
and ring particle to the initial amplitude of the ring particle's epicyclic motion. This 
gives the epicyclic amplitude of the centre of mass after one impact as 

s^^I^e.^Me, (4.22) 



mi + 7712 \nT'2 

Assuming that successive collisions are uncorrelated and occur stochastically with the 
mean time between consecutive encounters being r^e, the evolution of S is governed by 
the equation 









dt 


collisions 





(4.23) 



where = dN/dt can be expressed in terms of the surface density and an impact 
window function Wb+c{b) (see equation I4.12p . The quantity {£f) is the mean square 
value of Si for the ring particles. Using equation 14. 2[ this may be written in terms of the 
mean squares of the components of the velocity dispersion relative to the background 
shear, in the form 

(S'l) = Q-'iixl + 4{yi + 3Qxi/2Y). (4.24) 
We find in numerical simulations (ei) ~ 10^^ for almost all ring param eters. This value 



is ma inly determined by the coefficient of restitution (see e.g. figure 4 in lMorishima &: Salo 



20061). 



4.3.2 Stochastic excitation due to circulating particles 



Ring particles that are on circulating orbits, such as that given by path (a) in figure 
14. H exchange e nergy and angular momentum with the moonlet and therefore change 
its eccentricity. iGoldreich fc Tremaind ( 1l982l ) calculated the change in the eccentricity 
of a ring particle due to a satellite. We are interested in the change of the eccen- 
tricity of the moonlet induced by a ring of particles and therefore swap the satellite 
mass with the mass o f a rin g particle. Thus, rewriting their results (equation 64 in 
Goldreich fc Tremaine in our local notation without reference to the semi-major 

axis, we have 



AS', 



5.02 ml Q-^ 6" 



(4.25) 
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Supposing that the moonlet has very small eccentricity, it will receive stochastic im- 
pulses that cause its eccentricity to undergo a random walk that will result in in- 
creasing linearly with time, so that we may write 



dt 



/oo 
Waib)A£idil/U), (4.26) 
-oo 



circulating particles 

where d{l/tb) is the mean encounter rate with particles which have an impact parameter 
in the interval {b, b + db). Wa{b) is the window function describing the band of particles 
in circulating orbits. Setting d{l/tb) = 3T,Q\b\db / {2mi) , we obtain 

dS^ 



dt 



/oo 
Wa{b)m^\b\-^m^n-^db. (4.27) 
-oo 



circulating particles 

For high surface densities an additional effect can excite the eccentricity when gravita- 
tional wakes occur. The Toomre parameter Q is defined as 

where v is the velocity dispersion of the ring particles^. When the surface density is 
sufficiently high such that Q approaches unity, transient clumps will form in the rings. 
The typical leiigth s cale of these structures is given by the critical Toomre wave length 
fiDaisaka et"aDl2nnih 

At = Air^Gm-^, (4.29) 
which can be used to estimate a typical mass of the clump: 



rriT ~ ""[y) ^ = 47r^^'S=^^^^- (4-30) 

Whenever strong clumping occurs, mx should be used in the above calculation instead 
of the mass of a single particle mi. 



dSi 



dt 



/oo 
w;^{b)mT\b\-^j:G^n-^db, (4.31) 
-00 



circulating clumps 

where VF^(6) ~ Wa{b) is the appropriate window function. For typical parameters used 
in section 1121 this transition occurs at S ~ 200kg/m^. 



4.3.3 Equilibrium eccentricity 



Putting together the results from this and the previous section, we can estimate an equi- 
librium eccentricity of the moonlet. Let us assume the eccentricity 62, or the amplitude 
of the epicyclic motion S2, evolves under the influence of excitation and damping forces 
as follows 



dSi 
dt 



-2( 



7~ -\- T ~\- T 

e, collisions ' e, horseshoe ' e. 



+ 



del 



dt 



+ 



d£l 



collisions 



dt 



+ 



d£l 



dt ■ ^'-''^ 

circulating particles circulating clumps 

^The Toomre criterion used here was originally derived for a flat gaseous disc. To make use of it 
we replace the sound speed by the radial velocity dispersion of the ring particles. 
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The equilibrium eccentricity is then found by setting the above equation equal to zero 
and solving for £2. 

To make quantitative estimates we need to specify the window functions Waih), 
VI4+c(^) and Wd{h) that determine in which impact parameter bands the particles are 
in (see figure I4.ip . To compare our analytic estimates to the numerical simulations 
presented below, we measure the window functions numerically. Alternatively, one 
could simply use sharp cutoffs at some multiple of the Hill radius (see also section 
I4.6.ip . We already made use of this approximation as a simple estimate in the previous 
sections. The results may vary slightly, but not significantly. 

However, as the window functions are dimensionless, it is possible to obtain the di- 
mensional scaling of £2 by adopting the length scale applicable to the impact parameter 
to be the Hill radius th and simply assume that the window functions are of order unity. 
As already indicated above, all of the circularisation times scale as (xVLth G^^ S^^, 
or equivalently oc m2/(Sr|^fi). Assuming the ring particles have zero velocity disper- 
sion, the eccentricity excitation is then due to circulating clumps or particles and the 
scaling of d£l / dt due to this cause is given by equations 14.271 and 14.311 by 

^ oc m,rfT.G''n-\ (4.33) 

where rrii is either mi or m^. We may then find the scaling of the equilibrium value of 
£2 from consideration of equations 14.321 and 14.131 as 

£l oc "^rl. (4.34) 
m2 

This means that the expected kinetic energy in the non circular motion of the moonlet 
is oc mirjjQ'^ which can be viewed as stating that the non circular moonlet motion 
scales in equipartition with the mass moving with speed th^- This speed applies 
when the dispersion velocity associated with these masses is zero indicating that the 
shear across a Hill radius replaces the dispersion velocity in that limit. 

In the opposite limit when the dispersion velocity exceeds the shear across a Hill 
radius and the dominant source of eccentricity excitation is due to collisions, equation 
02] gives 

m2^2 = rni{£^) (4.35) 

so that the moonlet is in equipartition with the ring particles. Results for the two 
limiting cases can be combined to give an expression for the amplitude of the epicyclic 
motion excited in the moonlet of the form 

m2fl'^£^ = miQ^£^) + C,miQ^rjj, (4.36) 

where Cj is a constant of order unity. This indicates the transition between the shear 
dominated and the velocity dispersion dominated limits. 



4.4 Processes leading to a random walk in the 
semi-major axis of the moonlet 

We have established estimates for the equilibrium eccentricity of the moonlet in the 
previous section. Here, we estimate the random walk of the semi-major axis of the 



72 



moonlet. In contrast to the case of the eccentricity, there is no tendency for the semi- 
major axis to relax to any particular value, so that there are no damping terms and the 
deviation of the semi-major axis from its value at time t = grows on average as \/i 
for large t. 

Depending on the surface density, there are different effects that dominate the con- 
tributions to the random walk of the moonlet. For low surface densities, collisions and 
horseshoe orbits are most important. For high surface densities, the random walk is 
dominated by the stochastic gravitational force of circulating particles and clumps. In 
this section, we try to estimate the strength of each effect. 



4.4.1 Random walk due to collisions 

Let us assume, without loss of generality, that the guiding centre of the epicyclic motion 
of a ring particle is displaced from the orbit of the moonlet by Ai in the inertial frame, 
whereas in the co-rotating frame the moonlet is initially located at the origin with 
^2 = 0. When the impacting particle becomes bound to the moonlet, the guiding 
centre of the epicyclic motion of the centre of mass of the combined object, as viewed 
in the inertial frame, is then displaced from the original moonlet orbit by 

= - ^A. (4.37) 

mi + m2 m2 

This is the analogue of equation 14.221 for the evolution of the semi-major axis. For an 
ensemble of particles with impact parameter b, the average centre of epicyclic motion 
is (^i) = b. Thus we can write the evolution of A due to consecutive encounters as 



dA^ 



dt 



collisions 



-XiADr-,' (4.38) 

= ^ ri^^ 1^1' W,+c{b)db, (4.39) 

which should be compared to the corresponding expression for the eccentricity in equa- 
tion 1123 



4.4.2 Random walk due to stochastic forces from circulating 
particles and clumps 

Particles on a circulating trajectory (see path (a) in figure 14. ip that come close the 
moonlet will exert a stochastic gravitational force. The largest contributions occur for 
particles within a few Hill radii. Thus, we can crudely estimate the magnitude of the 
specific gravitational force (acceleration) due to a single particle as 



U = 7^- (4.40) 



Grrii 

When self-gravity is important, gravitational wakes (or clumps) have to be taken into 
account, as was done in section 14.3.21 Then a rough estimate of the largest specific 
gravitational force due to self gravitating clumps is 

{2rH + XtY 



fee ,^ , \ \2 ' (4-41; 
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where 2rH has been replaced by 2rH + A^- This allows for the fact that Ay could be 
significantly larger than rn, in which case the approximate distance of the clump to 
the moonlet is At (see also figure [3A] and discussion in section l3.2.ip . 

Following the formalism of chapter |3l we define a diffusion coefficient as D = 2tc{P), 
where / is an acceleration, t^. is the correlation time and the angle brackets denote a 
mean value. We here take the correlation time associated with these forces to be the 
orbital period, 27r^7^^. This is the natural dynamical timescale of the system and has 
been found to be a reasonable assumption from an analysis of the simulations described 
below. Using eqation 13. 4H we may estimate the random walk in A due to circulating 
particles and self-gravitating clumps to be given by 



dA^ 



dt 



dA^ 



AD^^n-^ = 1^7, n-^f I (4.42) 

2 

^2 



circulating particles 



16^ 7^) (4-43) 



dt 



circulating clumps 



AD,,^l-' = lQTi^l-''fi (4.44) 



2 



= IQi^n-'i A . (4.45) 

V(2r^, + AT)V 

This is only a crude estimate of the random walk undergone by the moonlet. In reality 
several additional effects might also play a role. For example, circulating particles 
and clumps are clearly correlated, the gravitational wakes have a large extent in the 
azimuthal direction and particles that spend a long time in the vicinity of the moonlet 
have more complex trajectories. Nevertheless, we find that the above estimates are 
correct up to a factor 2 for all the simulations that we performed (see below). 

4.4.3 Random walk due to particles in horseshoe orbits 

Finally, let us calculate the random walk induced by particles on horseshoe orbits. 
Particles undergoing horseshoe turns on opposite sides of the planet produce changes of 
opposite sign. Encounters with the moonlet are stochastic and therefore the semi-major 
axis will undergo a random walk. A single particle with impact parameter h will change 
the semi-major axis of the moonlet by 

= 2—^ — h. (4.46) 

mi + m2 

Analogous to the analysis in section 14.4. the time evolution of A is then governed by 
dA^ ' 



dt 



horseshoe 



6^ [ Sfi 161^ WJb)db. (4.48) 
mi J_oo 



Note that this equation is identical to equation 14.391 except a factor 4, as particles with 
impact parameter b will leave the vicinity of the moonlet at —b. 
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Name 




I 




r2 


dt 


Lx 


X Ly 






EQ5025 


50 


h 


r/m^ 


25 


in 


4s 


1000 


III 


X 


1000 


m 


7.2k 


EQ5050 


50 


h 


r/m^ 


50 


m 


4s 


1000 


III 


X 


1000 


m 


7.2k 


EQ20025 


200 


h 


l/ra^ 


25 


III 


4s 


1000 


III 


X 


1000 


m 


28.8k 


EQ20050 


200 


h 


l/ra^ 


50 


m 


4s 


1000 


111 


X 


1000 


m 


28.8k 


Eq40025 


400 


h 


r/m^ 


25 


m 


4s 


1000 


111 


X 


1000 


111 


57.6k 


EQ40050 


400 


h 


r/m^ 


50 


m 


4s 


1000 


m 


X 


1000 


m 


57.6k 


EQ40050DT 


400 


h 




50 


m 


40s 


1000 


III 


X 


1000 


m 


57.6k 


EQSOSODTW 


50 


h 


r/m^ 


50 


m 


40s 


2000 


III 


X 


2000 


m 


28.8k 


EQ20050DTW 


200 


h 


r/m^ 


50 


m 


40s 


2000 


III 


X 


2000 


m 


115.2k 


Eq40050DTW 


400 


h 




50 


m 


40s 


2000 


m 


X 


2000 


m 


230.0k 



Table 4.1: Initial simulation parameters. The second column gives the surface density 
of the ring. The third gives the moonlet radius. The fourth column gives the time 
step. The fifth and sixth columns give the lengths of the main box as measured in the 
xy-plane and the number of particles used, respectively. 



4.5 Numerical Simulations 

We perform realistic three dimensional simulations of Saturn's rings with an embedded 
moonlet and verify the analytic estimates presented above. The nomenclature and 
parameters used for the simulations are listed in table 14.11 



4.5.1 Methods 



We u se the GravTree code and a symplectic integrator for Hill's equations flQuinn et al. 
2010l ). Both are described in more detail in appendix |pl The gravitational forces are 
calculated with a Barnes-Hut tree flBarnes &: Hutlll986l ). 



To further speed up the calculations, we run two coupled simulations in parallel. 
The first adopts the main box which incorporates a moonlet. The second adopts an 
auxiliary box which initially has the same number of particles but which does not 
contain a moonlet. This is taken to be representative of the unperturbed background 
ring. We describe these so called pseudo shear periodic boundary conditions in more 
detail in appendix [Dl 

The moonlet is allowed to move freely in the main box. However, in order to prevent 
it leaving the computational domain, as soon as the moonlet has left the innermost part 
(defined as extending one eighth of the box size), all particles are shifted together with 
the box boundaries, such that the moonlet is returned to the centre of the box. This is 
possible because the shearing box approximation is invariant with respect to translations 
in the xy plane (see equation 14. ip . 

Collisions between particles are resolved using the i nstantaneous c o llision model and 
a velocity dependent coefficient of restitution given by [Bridges et al.l (119841 ): 



eiv 



mm 



0.34 



Icm/s 



-0.234 



(4.49) 



where v\\ is the impact speed projected on the axis between the two particles. 
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4.5.2 Initial conditions and tests 



Throughout this chapter, we use a distribution of particle sizes, ri, ranging from Im 
to 5m with a slope of g = —3. Thus dN/dri oc r^'^. The density of both the ring 
particles and the moonlet is taken to be 0.4g/c m^. This is at the low' er end of what 
has been assumed reasonable for Saturn's A ring (ILewis fc Stewartll2009[ ). The moonlet 
radius is taken to be either 50 m or 25 m. We found that using a larger ring particle 
and moonlet density only leads to more particles being bound to the moonlet. This 
effectively increases the mass of the moonlet (or equivalently its Hill radius) and can 
therefore be easily scaled to the formalism presented here. Simulating a gravitational 
aggregate of this kind is computationally very expensive, as many more collisions have 
to be resolved each time-step. 

The initial velocity dispersion of the particles is set to Imm/s for the x and y 
components and 0.4mm/s for the z component. The moonlet is placed at a semi-major 
axis of a = 130000 km, corresponding to an orbital period of P = 13.3 hours. Initially, 
the moonlet is placed in the centre of the shearing box on a circular orbit. 

The dimensions of each box, as viewed in the (x, y) plane, are specified in table W7[\ 
Because of the small dispersion velocities, the vertical motion is automatically strongly 
confined to the mid-plane. After a few orbits, the simulations reach an equilibrium 
state in which the velocity dispersion of particles does not change any more. 

We ran several tests to ensure that our results are converged. Simulation EQ40050DT 
uses a ten times larger time step than simulation EQ40050. Simulations EQ5050DTW, 
EQ20050DTW and EQ40050DTW use a box that is twice the size of that used in simulation 
EQ40050. Therefore, four times more particles have been used. No differences in the 
equilibrium state of the ring particles, nor in the equilibrium state of the moonlet have 
been observed in any of those test cases. 

A snapshot of the particle distribution in simulation EQ40050DTW is shown in fig- 
ure 1121 Snapshots of simulations EQ5025, EQ5050, EQ20025, EQ20050, EQ40025 and 
EQ40050, which use the smaller box size, are shown in figures 14. 3[ 14.41 and 14.51 In all 
these cases, the length of the wake that is created by the moonlet is much longer than 
the box size. Nevertheless, the pseudo shear periodic boundary conditions allow us to 
simulate the dynamical evolution of the moonlet accurately. 



4.6 Results 
4.6.1 Impact bands 

The impact band window functions Waih), Wh+dp) and Wdih) are necessary to estimate 
the damping timescales and the strength of the excitation. To measure those, each 
particle that enters the box is labelled with its impact parameter. The possible outcomes 
are plotted in figure 14.11 (a) being a circulating particle which leaves the box on the 
opposite side, (b) representing a collision with other ring particles close to the moonlet 
where the maximum distance from the centre of the moonlet has been taken to be twice 
the moonlet radius, (c) being a direct collision with the moonlet and finally (d) showing 
a horseshoe orbit in which the particle leaves the box on the same side that it entered. 

The impact band (b) is considered in addition to the impact band (c) because some 
ring particles will collide with other ring particles that are (temporarily) bound to the 
moonlet. Thus, these collisions take part in the transfer the energy and momentum 
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Figure 4.2: Particle distribution in the xy plane in simulation EQ40050DTW. The moon- 
let's wake is clearly visible. 



to the moonl et. Actually , if th e moonlet is simply a rubble pile of ring particles as 
suggested by iPorco et al.l (120071 ). then there might be no solid moonlet core and all 
collisions are in the impact band (b). 

Figures 14.61 and 14.71 show the impact bands, normalised to the Hill radius of the 
moonlet Th for simulations EQ5025, EQ20025, EQ40050, EQ5050, EQ20050 and EQ40050. 
For comparison, we also plot sharp cutoffs, at l.'bru and 2.5rH, as these have been used 
in same of the scaling laws presented above. 

Note that ^j^^ hcd^i ~ because ring particles are shifted from time to time when 
the moonlet is too far away from the origin (see section 14.5. ip . In this process, particles 
with the same initial impact parameter might become associated with the more than 
one impact band. 

Our results show no sharp discontinuity because of the velocity dispersion in the 
ring particles which is ~ 5mm/s. The velocity dispersion corresponds to an epicyclic 
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(a) S = 50kg/m2, ra = 25 m 




(b) E = 50kg/m2, = 50m 



Figure 4.3: Particle distribution in the xy plane in simulations EQ5025 and EQ5050 at 
a time when the moonlet and the ring particles are in equilibrium. The mean optical 
depth is r = 0.04. 
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(a) S = 200kg/m2, ra = 2.5 m 




(b) S = 200kg/m2, ra = 50 m 

Figure 4.4: Same as in figure but for simulations EQ20025 and EQ20050. The mean 
optical depth is r = 0.15. 
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(b) S = 400kg/m2, ra = 50 m 

Figure 4.5: Same as in figure but for simulations EQ40025 and EQ40050. The mean 
optical depth is r = 0.30. 
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amplitude of E^p ~ 20 m or equivalently ~ 0.4 r^^ and ~ 0.8 r^;/ (for r2 = 50m and 
r2 = 25m, respectively). Note that the Hill radius is approximately equal to the physical 
radius of the moonlet. This is because t he moonlet will always fill its Roche lobe with 
ring particles (see also lPorco et al.l 120071 ). 

The impact bands are almost independent of the surface density and depend only on 
the Hill radii and the mean epicyclic amplitude of the ring particles, or, more precisely, 
the ratio thereof. Thus we expect the impact bands to have a cutoff width of £rp- A 
simple approximation of the impact bands may then be given by 



Wa{h) 



( h/ru 



-2.5 



1 1 , 
— I — erf . 

2 2V ^vp/rH 

1 1 , 
— I — err 

2 2 



Wa{h) 



1.5 — h/ru 
-Wd{h). 



and 



(4.50) 

(4.51) 
(4.52) 



These curves are plotted on the right hand side of figures 14.61 and 14.71 being quali- 
tatively very similar to the measured impact bands. This approximation ignores the 
non-spherical shape of the Roche lobe and assumes = 
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Approximation 
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Figure 4.6: Impact bands for simulations with a moonlet radius r2 = 25 m. Left and 
right columns show the window functions measured in simulations and approximations 
of those, respectively (see equations l4.50tHl52|l . The surface density of the ring is, from 
top to bottom, E = 50kg/m2, E = 200kg/m2 and E = 400kg/m2. 
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Name 



Numerical results 



Analytic results 
collisions horseshoe 



EQ5025 

EQ5050 

EQ20025 

EQ20050 

EQ40025 

EQ40050 



18.9 
27.6 
9.0 
12.9 
2.8 
5.5 



18.4 
33.8 
3.9 
8.5 
2.1 
4.4 



1712 
3424 
428 
856 
214 
428 



Table 4.2: Eccentricity damping timescale Tg of the moonlet in units of the orbital 
period. The second column lists the simulation results. The third and fourth column 
list the analytic estimates of collisional and horseshoe damping timescales, respectively. 



4.6.2 Eccentricity damping timescale 

To measure the eccentricity damping timescale, we first let the ring particles and the 
moonlet reach an equilibrium and integrate them for 200 orbits. We then change the 
velocity of the moonlet and the ring particles within 2rH. The new velocity corresponds 
to an eccentricity of 6-10"'', which is well above the equilibrium value. We then measure 
the decay timescale Tg by fitting a function of the form 



The results are given in table 14.21 They are in good agreement with the estimated 
damping timescales from section 14.2.11 and I4.2.2[ showing clearly that the most impor- 
tant damping process, in every simulation considered here, is indeed through collisional 
damping as predicted by comparing equation 14.151 with equations 14.191 and 14.211 

4.6.3 Mean moonlet eccentricity 

The time averaged mean eccentricity of the moonlet^ is measured in all simulations after 
several orbits when the ring particles and the moonlet have reached an equilibrium state. 
To compare this value with the estimates from section H?T] we set equat ion 14 . 3 2 1 equal to 
zero and use the analytic damping timescale listed in table 14.21 The analytic estimates 
of the equilibrium eccentricity are calculated for each excitation mechanism separately 
to disentangle their effects. They are listed in the third, fourth and fifth column in 
table 14.31 The sixth column lists the analytic estimate for the mean eccentricity using 
the sum of all excitation mechanisms. 

For all simulations, the estimates are correct within a factor of about 2. For low 
surface densities, the excitation is dominated by individual particle collisions. For larger 
surface densities, it is dominated by the excitation due to circulating self gravitating 
clumps (gravitational wakes). The estimates and their trends are surprisingly accurate, 
as we have ignored several effects (see below). 

^Here, we are calculating the mean eccentricity, not the root-mean-square eccentricity. 



edit) = { 



eeg) + (6 ■ 10-^ - (ee,)) e-*/^^ 



(4.53) 



84 



Name 


Numerical results 


coiiioions 


Analytic 
circuidT;ing parxicieb 


results 

Circulating ciunips 


xoxai 


EQ5025 


6.1 • 10-^ 


4.5 • 10-« 


1.9 • 10-« 


0.3 • 10-« 


4.9 • 10-* 


EQ5050 


4.3 • 10-^ 


1.6 ■ 10~s 


1.4 • 10-« 


0.2 • 10-* 


2.1 • 10-* 


EQ20025 


6.4 • 10"^ 


4.6 • 10^* 


1.7-10-^ 


2.0 • 10"* 


5.3 • 10"* 


EQ20050 


4.9 ■ 10"^ 


1.6 ■ 10"^ 


1.4 ■ 10-* 


1.6 ■ 10"* 


2.7- 10"* 


EQ40025 


11.3 • 10"^ 


4.6-10-^ 


2.5 • 10"* 


8.4 • 10"* 


9.9 • 10"* 


EQ40050 


7.2 • 10-^ 


1.6- 10-« 


1.5 • 10-« 


4.9 • 10-* 


5.4 • 10-* 



Table 4.3: Equilibrium eccentricity (ceq) of the moonlet. The second column lists the simulation results. The third, fourth and fifth 
column list the analytic estimates of the equilibrium eccentricity assuming a single excitation mechanism. The last column lists the 
analytic estimate of the equilibrium eccentricity summing over all excitation mechanisms. 



Name 


Numerical results 


horseshoe 


collisions 


Analytic results 

circulating particles circulating clumps 


total 


EQ5025 


22.7m 


9.0m 


10.5m 


17.6m 


0.3m 


22.4m 


EQ5050 


12.6m 


4.5m 


4.8m 


4.4m 


0.1m 


7.9m 


EQ20025 


38.1m 


18.1m 


25.1m 


17.6m 


15.7m 


38.9m 


EQ20050 


22.7m 


13.6m 


9.6m 


4.4m 


4.8m 


14.7m 


EQ40025 


78.1m 


25.6m 


36.5m 


17.6m 


86.6m 


100.7m 


EQ40050 


45.6m 


12.8m 


13.4m 


4.4m 


31.4m 


36.7m 



Table 4.4: Change in semi-major axis of the moonlet after 100 orbits, A{At = 100 orbits). The second column lists the simulation 
results. The third till sixth columns list the analytic estimates of the change in semi-major axis assuming a single excitation mechanism. 
The last column lists the total estimated change in semi-major axis summing over all excitation mechanisms. 



4.6.4 Random walk in semi-major axis 

The random walk of the semi-major axis a (or equivalently the centre of epicychc motion 
A in the shearing sheet) of the moonlet in the numerical simulations are measured and 
compared to the analytic estimates presented in section U31 We run one simulation per 
parameter set. To get a statistically meaningful expression for the average random walk 
after a given time A{/S.t), we average all pairs of A{t) and A{t') for which t — t' = At. 
In other words, we assume the system satisfies the Ergodic hypothesis. 

^(At) then grows like y/At and we can fit a simple square root function (see chapter 
|3]). This allows us to accurately measure the average growth in A after At = 100 orbits 
by running one simulation for a long time and averaging over time, rather than running 
many simulations and performing an ensemble average. The measured values are given 
in the second column of table l4.4l The values that correspond to the analytic expressions 
in equations 14. 48^ 14. 39^ 14.431 and 14.451 are listed in columns three, four, five and six, 
respectively. 

For low surface densities, the evolution of the random walk is dominated by collisions 
and the effect of particles on horseshoe orbits. For large surface densities, the main effect 
comes from the stochastic gravitational force due to circulating clumps (gravitational 
wakes) . 

The change in semi-major axis is relatively small, a few tens of metres after 100 
orbits (= 50 days). The change in longitude AA that corresponds to this is, however, 
is much larger. 

Using the results from chapter [3], namely equations 13.561 and I3.60[ we can predict 
the long term evolution of A. In the stochastic and laminar case, (AA)^ grows like ~ t^ 
and ~ t^, respectively. This behaviour allows us to discriminate between a stochastic 
and laminar migration by observing AA over an extended period of time. 

Results from individual simulations (i.e. not an ensemble average) are plotted in 
figure Sin Here, A A is expressed in terms of the azimuthal offset relative to a Keplerian 
orbit. One can see that for an individual moonlet, the shift in azimuth can appear to 
be linear, constant or oscillating on short timescales (see curve for simulation EQ20025). 
This is partly because of the lower order terms in equation I3.55[ though, on average, 
the azimuthal offset grows very rapidly, as ~ t^^^ for t ^ r (see equation 13.561) . 

4.7 Conclusions 

In this chapter, we have discussed the dynamical response of an embedded moonlet in 
Saturn's rings to interactions with ring particles both analytically and by the use of 
realistic three dimensional many particle simulations. Both the moonlet and the ring 
particle density were taken to be 0.4g/cm^. Moonlets of radius 25 m and 50 m were 
considered. Particle sizes ranging between 1 m and 5 m were adopted. 

We estimated the eccentricity damping timescale of the moonlet due to collisions 
with ring particles and due to the response of ring particles to gravitational perturba- 
tions by the moonlet analytically. We found the effects due to the response of particles 
on horseshoe and circulating orbits are negligible. On the other hand the stochastic 
impulses applied to the moonlet by circulating particles were found to cause the square 
of the eccentricity to grow linearly with time as did collisions with particles with non 
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Figure 4.8: Offset in tfie azimutfial distance due to tlie random walk in tfie semi- 
major axis a measured in simulations EQ5025, EQ5050, ECj20025, EQ20050, EQ40025 
and EQ40050. Individual moonlets may show linear, constant or oscillatory growth. On 
average, the azimuthal offset grows like t^^'^. 



zero velocity dispersion. A balance between excitation and damping processes then 
leads to an equilibrium moonlet eccentricity. 

We also estimated the magnitude of the random walk in the semi-major axis of 
the moonlet induced by collisions with individual ring particles and the gravitational 
interaction with particles and gravitational wakes. There is no tendency for the semi- 
major axis to relax to any particular value, so that there are no damping terms. The 
deviation of the semi-major axis from its value at time t — grows on average as y/i 
for large t. 

From our simulations we find that the evolution of the eccentricity is indeed domi- 
nated by collisions with ring particles. For large surface densities (more than 200kg/m^) 
the effect of gravitational wakes also becomes important, leading to an increase in the 
mean steady state eccentricity of the moonlet. When the particle velocity dispersion 
is large compared to ^th, we obtain approximate energy equipartition between the 
moonlet and ring particles as far as epicyclic motion is concerned. 

Similarly, the random walk in the semi-major axis was found to be dominated by 
collisions for low surface densities and by gravitational wakes for large surface densi- 
ties. In addition we have shown that on average, the difference in longitude A A of a 
stochastically forced moonlet grows with time like t^^"^ for large t. 

The radial distance travelled within 100 orbits (50 days at a distance of 130000km) 
is, depending on the precise parameters, of the order of 10-lOOm. This translates to 
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a shift in longitude of several hundred kilometres. The shift in A is not necessarily 
monotonic on short timescales (see figure H78|) . We expe ct that such a s h ift sh ould be 
easily detectable by the Cassini spacecraft. And indeed, iTiscareno et al.l (120101 ) report 
such an observation, although the data is not publicly available yet, at the time when 
this thesis was submitted. 

There are several effects that have not been included in this study. The analytic 
discussions in sections 14. 2[ 14.31 and 14.41 do not model the motion of particles correctly 
when their trajectories take them very close to the moonlet. Material that is tem- 
porarily or permanently bound to the moonlet is ignored. The density of both the ring 
particles and the moon let are at the lower end of what is assumed to be reasonable 
fiLewis fc Steward l2009h . which has the advantage that for our computational setup 
only a few particles (with a total mass of less than 10% of the moonlet) are bound to 
the moonlet at any given time. 
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Chapter 



Planetesimal formation in a turbulent 

disc 

To explain all nature is too difficult a task for any 
one man or even for any one age. 'Tis much better to 
do a little with certainty, & leave the rest for others 
that come after you, than to explain all things by 
conjecture without making sure of any thing. 

Isaac Newton, unpublished preface to Opticks, 1704 

In chapter [31 we discussed the evolution of planets and planetary systems under the 
effect of stochastic forces resulting from turbulence. The stochastic migration of moon- 
lets in Saturn's rings was discussed in chapter HJ In this chapter, we look at another 
system, in which turbulence is important, namely the formation of planetesimals, the 
building blocks of planets. 

The formation mechanism of planetesimals in proto-planetary discs is hotly debated. 
Currently, the favoured model involves the accumulation of metre-sized objects within 
a turbulent disc. These high spatial densities then promote the formation of planetes- 
imals, for example via gravitational instability as discussed in more detail in chapter 
[U In numerical simulations that model this process, one can at most simulate a few 
million particles as opposed to the several trillion metre-sized particles expected in a 
real proto-planetary disc. Therefore, single particles are often used as super-particles 
to represent a distribution of many smaller particles. It is assumed that small-scale 
phenomena do not play a role and particle collisions are not modelled. However, the 
super-particle approximation can only be valid in a collision-less or strongly collisional 
system. In many recent numerical simulations this is not the case. 

Here, a scaled system is studied that does not require the use of super-particles. We 
compare this new method to the standard super-particle approach. We find that the 
super-particle approach produces unrehable results that depend on numerical artifacts 
such as the gravitational softening in both the requirement for gravitational collapse and 
the resulting clump statistics. Our results show that short-range interactions (collisions) 
have to be modelled properly. 

The outline of this chapter is as follows. In section 15. we review the important 
timescales in the problem and show that the super-particle approach breaks down in 
high density regions. In section 15. 2[ we discuss both the super-particle and the scaled 
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particle method. The numerical implementation and the initial conditions of the sim- 
ulations are presented in section 15.31 Results of numerical simulations using both the 
super-particle approach and our scaling method are presented in section 15.41 We con- 
clude with presenting resolution constraints for numerical simulations and discuss the 
implications for the planetesimal formation process through gravitational instability in 
section 15. 5[ 



5.1 Orders of magnitude 
5.1.1 Definition of timescales 

The aim of this work is to simulate the dynamics of dust (or particles) interacting 
gravitationally inside a turbulent proto-planetary disc. Each particle is subject to five 
different physical processes, each operating on a typical timescale: 

Stopping time Tg 

Each particle of size a feels the effects of the surrounding gas through a linear 
drag force. This force can be written as Fdrag = — Vp), where is the 

stopping time of a parti cle of mass rup, Vn is the velocity of the gas, and Vp is the 



the degree of coupling to the gas (e.g. IChiang fc Youdinll2009[ ) 



velocity of the particle f lWeidenschillind Il977al). The stoppiii g time depends on 




a < 7A Epstein drag 

Tg = J Pa"" 2 ~ (5.1) 

4ppa Stokes drag. 



where pp and pg are the internal density of the particles and the density of the 
gas, respectively. A is the mean free path of gas particles and the sound speed. 

Physical collision timescale Tc 

Dust particles will suffer a physical collision with another dust particle on a 
timescale of Tc = {acVpn)~^, where cTc is the geometrical cross-section of the parti- 
cles, Vp is the particle velocity dispersion, and n is the number density of particles 
in a homogeneous medium. 

Orbital timescale Tg 

The timescale associated with the particles orbiting the central object in a local 
shearing patch is the epicyclic period Tg = 2-71^^^, where fl is the angular velocity 
at the semi-major axis of the shearing box. 

There are two limits for the gravitational interactions between dust particles, long and 
short range. The long-range interaction can be seen as a collective process involving 
interactions between clouds of particles. On the other hand, the short-range interaction 
is important for resolving close approaches between pairs of particles, i.e. gravitational 
scattering. For the sake of clarity, we separate these two processes. 

Gravitational collapse timescale tqi 

On a length scale A ^ 5.,.-, where 5r is the average distance between neighbouring 
dust particles, the long-range interaction of the dust particle distribution can be 
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approximated by a continuous density p. The gravitational collapse timescale, 
which is of the order of the free fall time of the system, can be defined by 
Tgi = 1/ y/Gp. 



Gravitational scattering timescale tgs 

Short-range gravitational interactions are interactions between a pair of single 
particles that can result in a scattering event. As with physical collisions, the 
timescale for a gravitational interaction depends on the cross-section, tgs = 
(o"GVp?T.)~^. Here, ac ~ G'^nip/Vp is the gravitational cross-section of the particle. 
Note that the cross-section is velocity dependent, which makes it qualitatively 
different from physical (billard ball) collisions between two particles. 

Particles might also be affected by excitation from a turbulent background state on 
a timescale defined by the turbulence itself. Unlike in chapter [3l we here model the 
excitation in an even more simplified way, in which the turbulence is scale independent 
and has no timescale associated with it. 



5.1.2 The real physical system 



To identify the dominant physical processes in a proto-planetary disc, one has to quan- 
tify and compare the relevant timescales as defined above. In the following discus- 
sion and in the numerical simulations, we assume Rq = 1 AU, a gas disc thickness of 
H/Rq = 0.01, and a minimum mass solar nebula (MMSN) with surface gas density^ 
E = 890 g cm~^. One usually assumes a solid to gas ratio of ps/ Pg = 0.01. We, however, 
assume a solid to gas ratio of unity. This value is justified by numerical simulations 
( jJohansen et al.l 120071 . |2009[ ) that demonstrate over- densities of the order of 10^ — 10^ 
can easily occur because of the interaction with a turbulent gas disc, the streaming 
instability, or vertical settling. Thus, we are only interested in the late stages of the 
planetesimal formation process. The details of how the initial over-density has been 
achieved are not relevant for the discussion presented in this chapter. We also note 



that significantly different models of the solar nebula have been proposed ( lDesclj|2007| ). 
However, all results can easily be scaled to different scenarios. 

We simplify the system by assuming that all solid components of the disc are in 
metre-sized particles. We assume that these boulders have a typical velocity dispersion 
caused by gas turbulence v^ 



O.OS c.s ^ 30m/s, where c ., is the local sound speed, in 



accordance with numerical results (iJohansen et al.l 120071 ). Since Vp <^ Cg, the particles 



tend to sediment toward the mid-plane, forming a finite thickness dust layer due to the 
non-zero velocity dispersion. 

For metre-sized boulders, the physical cross-section is cTc ~ Im^, whereas the grav- 
itational scattering cross-section is ac ~ 10^^^ m^ showing clearly that gravitational 
scattering is irrelevant to the dynamics of dust particles embedded in a disc. However, 
all other timescales are roughly equivalent. Metre-sized boulders are weakly coupled to 
the gas with a stopping time ~ Te ~ fl~^ (IWeidenschillingI Il977al ) . The physical col- 
lision timescale is Tc = 7.3 and the long-range gravitational interaction timescales 
is Tgi ~ 



^Note that there are different values for MMSN in the hterature. 
smaller than the one used in chapter [S] 



This value here is a factor 5 
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This physical syste m is expected to become gravitationally unstable, according to 
the Toomre criterion (IToomrd Il964| ). The instability occurs when the gravitational 
collapse timescale tqi is shorter than the transit time due to random particle motions 
A/vp and the orbital timescale. Assuming that the particles can be modelled by an 
isotropic gas^ with a sound speed Vp, the system becomes unstable when Q ~ 1 (equa- 
tion |12H])- In that case, the most unstable wavelength is given by equation I4.29[ In 
the following, we compare the physical parameters from this section to their counter- 
part in numerical simulations. We first summarise the super-particle approach before 
presenting the scaling method. 



5.2 Methods 

Two different kinds of simulation are considered in this chapter: 

1. Particles are assumed to be point masses and have no physical size, the gravita- 
tional field of the particles being approximated with a smoothing length to avoid 
numerical divergences (see section 15.3.11) . 

2. Particles have a physical size and, therefore, no gravitational smoothing is required 
but physical collisions must be included. 

We refer to the particles as super-particles and scaled particles, respectively. 



5.2.1 Super-particle approximation 

Let us quickly review the standard approach using super-particles before moving on to 
studying our favoured method that uses scaled particles. 

A super-particle represents many smaller particles. The dynamics of the small 
particles are not calculated exactly. It is assumed that all particles behave similarly and 
in a collective manner. To simulate gravitational interactions between super-particles, 
one has to use a softened potential. Without this, the gravitational scattering cross- 
section of super-particles becomes too large and super-particles undergo gravitational 
scattering events, which are unphysical since these events never occur in a real system. 
Individual particle-particle collisions are not modelled. 

There are various examples where this approach is used successfully. For example, 
smoothed particle hydrod ynamics (SP H) uses the super-particle approximation to sim- 



ulate many gas molecules (jLucylll977l ). These systems are often assumed to be strongly 
collisional to ensure thermodynamical equilibrium inside each super-particle. One can 
therefore assign collective properties to clouds of particles such as pressure and tem- 
perature. Another example is the evolution of galaxies. When two galaxies collide, 
individual particles (stars) will usually not undergo gravitational scattering events or 
physical collisions. In that case, the super-particle approach models a collision-less 
system in which collective dynamics are the only important physical process. In both 
cases, the super-particle approximation is a valid approach for simulating the system 
numerically, but will break down as soon as the system is marginally collisional. 

As an example, we consider two clouds of particles undergoing a "collision" . In the 
strongly collisional case, the clouds will slowly merge and the thermodynamic variables 



^This is formally not valid in the studied regime, as it is not strongly collisional. 
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(e.g., temperature) will diffuse between the clouds, the particles inside each cloud fol- 
lowing a random walk trajectory due to numerous collisions. On the other hand, in 
the collision- less regime, the two clouds simply do not see each other because there is 
no short-range interaction present between the particles. In a marginally collisional 
regime, some particles will collide with particles of the other cloud, leading to a partial 
thermalisation of the velocity distribution, but some other particles will not have any 
collision at all and will follow approximately a straight line. Evidently, the outcome of 
this event cannot be described using a super-particle approach. 



5.2.2 Scaling method 

The idea of the scahng method is that one should keep all important timescales in a 
numerical simulation as close to those of the real physical system as possible and model 
all particle collisions explicitly. 

The numerical system consists of iVnum particles in a box with length if, simulating 
a patch of the disc at a fixed radius Rq. The particle mass is mnum = ^H'^/Nnum {H 
is the box size and one scale height). Thus, the density in the box and the long-range 
gravitational interactions are unchanged compared to the real system. 

The gravitational scattering cross-section of the numerical system is then given by 

num p 

For the initial surface density and velocity dispersion used in the simulation (see above), 
one finds 



0.0001 AU^ 



Owl 



N2 _/V2 
num num 



(5.3) 



The physical collision cross-section in the simulation is crc,num = '^O'nxim where Onum, 
is the radius of the particles in the simulation. We derive two constraints from the 
physical and gravitational collision cross-sections: 

1. Same mean free path in simulation and real physical systems 

That means A''numO'c,num = Nac, where cxc and N are the geometrical cross-section 
and the number of particles in a region of size in a real disc, respectively. 
In other words, this condition ensures that the physical collision timescale in the 
simulation is exactly the same as in the real system. 

2. Negligible gravitational scattering 

When two particles approach each other, the outcome should be a physical colli- 
sion, which means that (7c,num ^ co.num- The gravitational scattering timescale 
remains long compared to the physical collision timescale. 

The first condition places a constraint on the particle size in the simulation. With 
N — 4.7 ■ 10^* and a particle radius of a = 1 m as found in a real disc within a box of 
volume H^, using the parameters from above one finds that 



^ = 4£11aU. (5.4) 

TT V A^numTT V^num 
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We note that once this condition is satisfied in the initial conditions it will be automat- 
ically satisfied at all times. In simulations, we typically find that the collision timescale 
is reduced by more than one order of magnitude during the gravitational collapse. 

The second condition is then satisfied by changing the number of particles. An 
interesting result is that if iVnum > 10, the second condition is easily satisfied if the 
first one is satisfied. However, acnum depends strongly on the velocity dispersion. In 
particular, a velocity dispersion 5 times smaller than the initial value (as found in some 
simulations) leads to an increase by a factor of 600 in the gravitational scattering cross- 
section. Therefore, we suggest using a large safety factor (A^num > 10^) to ensure that 
the second condition is always satisfied, even for significantly smaller Vp. This condition 
also allows us to estimate when our approach breaks down, namely when the number 
of clumps in the system is so low (A^ < 10 ~ 100) that gravitational scattering becomes 
important. 

Although it would be helpful as a further simplification, it is not possible to perform 
the above mentioned simulation in two dimensions. In that case, the filling factor, which 
is defined to be the ratio of the volume (or area) of all particles to the total volume (or 
area), is of the order of one for the above parameters. We note that in 2D an increase 
in particle number (and decrease in particle size as required by equation 15. 4p does not 
decrease the filling factor if the collisional lifetime remains constant. 

5.3 Numerical simulations 

We perform our simulations in a cubic box with shear periodic boundary conditions in 
the radial (x) and periodic boundary conditions in the azimuthal (y) and vertical (z) 
directions, as described in appendix |Dl In the local approximation, the force per mass 
on each particle is a sum of the contributions from Hill's equations (see equations ID.81 - 
ID.lOp and interaction terms. The interaction terms -Fint are divided into components 
related to self gravity (including vertical gravity), physical collisions between particles, 
or drag and excitation forces: 

-^int -^grav ~l~ -^col ~l~ -^drag ~l~ -^turb- (^■^) 

We solve the resulting equations of motion with a leap frog (kick drift kick) time- 
stepping scheme. In the following subsections, we describe how each of these terms are 
calculated and what their physical relevance is. 

The box size of 0.01 AU was chosen such that there are always several unstable modes 
At in the box, when the system is pushed into the regime of gravitational collapse, as 
estimated by equation 14.291 

5.3.1 Self gravity 

In the A^-body problem that we consider, we have to solve Newton's equations of uni- 
versal gravitation for a large number of particles A^. Calculating the gravity for each 
particle from each other particle results in O {N'^) operations. To reduce the number 
of operations, we use an approximation, namely a Barnes-Hut (BH) tree code (see 
appendix [D|) . 

We use 8 rings of ghost boxes in the radial and azimuthal direction. The system 
that we are interested in is only marginally gravitationally stable and due to the finite 
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number of ghost boxes, we introduce a slight asymmetry that will become important 
after several orbits. In the beginning of our simulations, we integrate for many orbits 
in order to reach a stable equilibrium. We then push the system into a gravitationally 
unstable regime (see section l5.3.5l for details). During the stable phase, horizontal over- 
concentrations are to be expected and are indeed observed because of this asymmetry. 
However, since in a real proto-planetary disc the gravitational instability is considered 
to appear locally (e.g., inside a vortex or a similar structure), this effect of having a 
preferred location for gravitational instability is not unphysical and does not affect any 
conclusions made. We tested this by applying a linear cut-off to the gravitational force 
at a distance of one box length (A. Toomre, private communication). In this case, there 
was no preferred location in the box and the simulation evolved in exactly the same 
way. 



5.3.2 Physical collisions 

Physical collisions between particles are treated by checking whether any two particles 
are overlapping after each time-step. Using the already existing tree structure from 
the gravity calculation, one can again reduce the computational costs from O {N"^) to 
O (NlogN). There is a small probability that we might miss a collision if the time-step 
is too large. To ensure that our results are converged, we tested that a reduced time- 
step does not alter the outcome. See appendix [D] for more details on how collisions are 
resolved. 

The particle radius is given by equation 15.41 In order to keep the collision timescale 
Tc close to unity, which is numerically the worst case scenario because all timescales are 
equivalent, we increase the radius by a factor of 4 in all simulations. 



5.3.3 Drag force 

Each particle feels a drag force. The background velocity of the gas Vg is assumed to 
be a steady Keplerian profile 

3 

Vg = -^^x ey. (5.6) 

Although accretion flows are turbulent, we choose to use a simple velocity profile so 
that the behaviour of self-gravitating particles can be understood in conditions that 
can be easily controlled. 



5.3.4 Random excitation 



The system of particles described above is dynamically unstable even without self- 
gravity, as the coefficient of restitution and the stopping time do not depend on the 
particle velocities. The particles can either settle down in the mid-plane and create a 
razor-thin disc if the stopping time is too short, or they can expand vertically forever if 
the stopping time is above some critical value and the excitation mechanism is provided 
only by collisions. Dust particles in accretion discs are found in the settling regime. 
However, a complete settling never occurs in accretion discs b ecause the backgroun d 
flow is always turbulent due to the Kelvin- Helmholtz i nstabilitv (IJohansen et al.ll2006l) . 
the MRI, or other hydrodynamic instabilities (see e.g. iLesur fc Papaloizod 120091 ) which 
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diffuse particles vertically (iFromang &: Papaloizoul 120061 ) . This stirring process of tur- 



bulence is not present in the gas velocity field of the simulations presented here (see 
equation 15.61) . To approximate the turbulent mixing, we add a random excitation 
(white noise in space and time) to the particles in the simulation. This allows us to 
have a well defined equilibrium in which the system is stable rather than starting from 
unstable initial conditions that might influence the final state. 

We perturb the velocity components of each particle after each time-step by Avj = 
\/Si^, where 6t is the current time-step and ^ is a random variable with a normal distri- 
bution around and variance s. This excitation mechanism heats up the particles and 
allows us to have a well defined three dimensional equilibrium as shown in section 15.41 
In our simulations, we use a value of s = 1.3 ■ 10^^ m^s~^ and s = 8.9 ■ 10^^ m^s~^ for the 
simulations without and with collisions, respectively. These values were chosen such 
that the equilibrium state is approximately equivalent for both types of simulations. 

In comparison to chapter [3l here we try to simplify things as far as possible. Thus, 
the description of stochastic forces used leads to uncorrelated and scale independent 
noise. 



5.3.5 Initial conditions 

All particles, which have equal mass, smoothing length, and physical size, are placed 
randomly inside the box in the xy plane. We note that particles have either a physical 
size or a smoothing length associated with them, depending on the type of simulation 
we perform. In the z-direction, the particles are placed in a layer with an initially 
Gaussian distribution about the mid-plane and a standard deviation of 0.05 H. We 
allow the system to reach equilibrium by integrating it until t = 30^2"^ (4.8yrs). 

We tested various ways of pushing the system into the unstable regime. We here 
present two different scenario. First, one can switch off the excitation mechanism. The 
systems starts to cool down and thus becomes unstable. Second, one can shorten the 
stopping time, which leads to enhanced cooling. 

Following one of the above descriptions, the system becomes gravitationally unstable 
and bound clumps form within a few orbits. We did not see any qualitative differences 
as long as we started from a well defined equilibrium state. This shows that, even 
though we simplified the treatment of turbulent stirring and gas drag, the results are 
generic. 



5.4 Results 

5.4.1 Super particles 

We first present simulations without physical collisions (super-particles) which rely on 
the smoothing length b to avoid any divergences. Although we use a tree code and 
a smoothed potential for the force calculation is assumed, the results are equivalent 
to an FFT-based method (Fast Fourier Transform) where the grid length acts as an 
effective smoothing length. In general, to check the numerical resolution, the particle 
number A^num is increased and the smoothing length b is reduced independently. A 
simulation is resolved when the result is independent of both A'num and b. This turns 
out to be impossible in the present situation. The main reason is that the smallest 



96 



scale in a gravitational collapse will ultimately depend on the smoothing length. By 
varying both parameters at the same time, an empirical scaling of 6 ~ 1/ V^num works 
fairly well if A^num is large enough. The square root in this scaling comes from the fact 
that we are basically treating a two dimensional system. However, this procedure is not 
justified and is unphysical. The smoothing length was introduced to avoid divergent 
terms and not to model any small-scale physical process. Therefore, it should not have 
any impact on the physical outcome of the simulation. Incidentally, the existence of 
this dependency indicates that short-range interactions are important for the result of 
these simulations and should be modelled with care. 



In figure 5.1(a) , we plot the velocity dispersion as a function of time. The simulations 
begin from the stable equilibrium described in section 15.21 After t = 30 we 
switch off the excitation mechanism. The system continues to cool and it becomes 
gravitationally unstable within one orbit, as estimated by equation 14. 281 Once clumping 
occurs, the velocity dispersion begins to rise again. All simulations use the particle 
radius given by equation 15.41 as a smoothing length except the ones labelled LS, which 
use a value ten times larger. The larger softening length is approximately 128-th of the 
box length and illustrates the kind of evolution expected from an FFT method using a 
128'^ grid. This is the resolution used in other studies, such as I Johansen et al.l ( 120071 ). 



The velocity dispersion evolution is also plotted in figure 5.2(a) for those simula- 



tions, in which the damping timescale has been decreased while keeping the excitation 



unchanged. The results are qualitatively very similar to those presented in figure 5.1(a) 



Snapshots of the particle distribution of two simulations are shown in figure 15. 3[ 
Both simulations use 160 000 particles and all parameters, except the smoothing length, 
are the same. The top row is the simulation with a smoothing length given by equation 
15.41 whereas the bottom row uses a smoothing length that is ten times larger. The 
simulations correspond to the blue (medium dashed) and red (solid) curve in figure 



5.1(a) We show the snapshots to illustrate the importance of the smoothing length in 
simulations without physical collisions. One can see that the simulations differ already 
before clumps form. Stripy structures appear on a scale given by equation 14.291 As 
soon as we enter the unstable regime at t ~ 32^"^ (i.e., the velocity dispersion begins 



to rise, see also figure 5.1(a) ) the simulations evolve very differently. The simulation on 
the top row of figure ESI forms many clumps at an early time, whereas the bottom row 
simulation forms only a few, more massive clumps at later times. 

This can be confirmed by looking at spectra of the same two simulations as shown in 



figure 5.5(a), These spectra were generated by mapping the particles onto a 128 x 128 



grid in the xy plane and computing the fast Fourier transform of the map ping in the x 



direc tion. The resulting spectra were finally averaged in the y direction (see lTanga et al. 



2004 for a complete description of the procedure). As suggested by the snapshots, the 
nonlinear dynamics of the gravitational instability strongly depends on the smoothing 
length used. In particular, one observes that smaller scales are amplified more slowly for 
larger smoothing lengths. The resulting spectra at t = 40Q~^ also differ significantly. 
With a small enough smoothing length, the spectrum looks almost flat, whereas a large 
smoothing length introduces a cutoff at k/2'ir ~ 10. We also note that the smoothing 
length in the latter case is of the order of the grid size (/c/27r ~ 100). The cutoff observed 
clearly demonstrates that the smoothing length modifies the dynamics on scales up to 
10 times larger than the smoothing length itself. 
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(a) Super-particles. The simulation labelled LS has a ten times larger smoothing length. 
The curves do not overlap because the simulations have not converged. 
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(b) Scaled-particles. All curves overlap at early times because the simulations are 
converged (see text). 

Figure 5.1: The velocity dispersion as a function of time in different runs as an indicator 
of convergence. Simulations use super-particles (top) and scaled-particles (bottom). 
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(a) Super-particles. The simulation labelled LS has a ten times larger smoothing length. 
The curves do not overlap because the simulations have not converged. 
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(b) Scaled-particles. All curves overlap because the simulations are converged. 



Figure 5.2: Same as in figure \57[\ but here the damping timescale is decreased while the 
excitation is unchanged. This leads to qualitatively very similar results. 
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Figure 5.3: Super-particles: Snapshots of the particle distribution in the xy plane. Both simulations use 160 000 particles with different 
smoothing lengths. The simulation on the bottom uses smoothing length ten times larger than the one on the top. The snapshots were 
taken (from left to right) at t = 0, 30, 37, 40 With a large smoothing length, the outcome looks very different, the system is more 
stable, more stripy structure can be seen and clumps form later, if at all. 



o 




Figure 5.4: Scaled particles: Snapshots of the particle distribution in the xy plane. The simulations (from top to bottom) use 40 000, 
320 000, 640 000 particles with their physical size given by equation 15.41 In all simulations, the collision timescale was kept constant. 
The snapshots were taken (from left to right) at t = 0,30,37,40^2"^. The simulation with 40 000 particles has a large filling factor by 
the last frame, which prevents clumps from forming. The intermediate resolution simulation (middle row) and the highest resolution 
simulation (bottom row) have more and smaller particles and thus a smaller filling factor. The results (i.e., number of clumps in the 
last frame) are very similar in the intermediate and high resolution simulations. 



5.4.2 Scaled particles 



The velocity dispersion evolution of simulations with 

particles including physical collisions is presented in figure 5.1(b) and figure 5.2(b) 



160 000,320 000,640 000 

In 



these runs, no gravitational smoothing length is needed. Again, the simulations start 
from a stable equilibrium and after t = 30 we change the same parameters as in 
the non-collisional case to push the system into the gravitationally unstable regime. In 
we increase the damping timescale, whereas we switch off the excitation in 
Snapshots of the particle distribution for three runs are plotted in figure 



figure |5.1(b) 



figure |5.2(b) 



The middle and bottom rows of snapshots correspond to the purple (medium 



dashed) and dark blue (long dashed) lines of figure 5.1(b) 



One can see in figures 5.1(b), 5.2(b) as well as figure l5^ that the intermediate and 



high resolution runs produce very similar results. We call those simulations converged, 
as a change in particle number does not change the outcome. An additional, lower 
resolution run (A^ = 40 000) has not converged because the filling factor is too large in 
dense regions, preventing clumps from being gravitationally bound. The problem does 
not occur in the non-collisional cases because particles are allowed to overlap. Since 
the filling factor scales with particle number as 1/ A/A^num, this issue is resolved for runs 
with more than a few hundred thousand particles. We note that the velocity dispersion 
might differ at later stages for the converged runs because the final clump radius is still 
determined by the physical particle radius and therefore by the particle number. This 



can be seen in particular in figure 5.1(b) where the velocity dispersion of the 160 000 
run reaches an equilibrium value before the runs with larger particle number. Note 
however, that up to that point, the results are converged. 



We show the density distribution spectrum in figure 5.5(b) as a function of the 
number of particles. The resulting spectra do not depend on the number of particles 
in the system, as expected. In comparison with figure 5.5(a), the simulations with 



collisions are converged since the dynamic has not controlled by the numerical parameter 
A^, which is the only free parameter in the system. 

At very late times (t > 38^2^^), the spectra begin to show a systematic trend towards 
more structure on smaller scales for runs with larger A^. This is due to the filling factor, 
already mentioned above. When the size of a clump in the simulation is only determined 
by the size of the individual scaled particles that it consists of, our approach breaks 
down. One can also see these clumps forming by looking at the velocity dispersion in 
figure 5.1(b) The particles inside a clump begin to dominate the velocity dispersion 
over the background after t ~ 38fi~^. A clear indication of this is the spiky structure 
with a typical correlation time of ~ corresponding to different clumps interacting 
and merging with each other. One way ar ound this issue, which w ill be considered in 



future work, is to allow particles to merge flMichikoshi et al.ll2009l ). Using this specific 
accretion model, the mass of the clump can then be used as an upper limit of the 
clump mass. Eventually, one should also worry about possible effects leading to the 
destruction of an already formed clump. 
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(a) Super-particles: Spectrum with 160 000 particles and two different smoothing 
lengths that differ by a factor 10 at times t = 30 (red, solid curve), t — 36 (green, 
long dashed curve), and t = 38 (blue, short dashed curved). The spectra at similar 
times are not on top of each other because the simulations have not converged. 
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(b) Scaled particles: Sp ectrum with 160 000, 320 000, and 640 000 particles (colours 
are equivalent to figure 5.5(a)). At each time, the spectra are converged, i.e., the 
spectra are on top of each other and independent of the particle number, besides 
the noise level on very small scales. At late times {t > 38, blue curve), one begins 
to see more structure on small scales in simulations with larger number of particles 
(see text). 



Figure 5.5: Density distribution spectrum 
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5.5 Discussion and implications 



In this chapter, we have shown that convergence in A^-body simulations of planetesimal 
formation via gravitational instability can be achieved when taking into account all rel- 
evant physical processes. It is absolutely vital to simulate gravity, damping, excitation, 
and physical collisions simultaneously, as the related timescales are of the same order 
and therefore all effects are strongly coupled. 

A set of simulations is defined to be converged when the results do not depend on the 
particle number or any other artificial numerical parameter such as a smoothing length. 
As a test case, a box with shear periodic boundary conditions was used, containing 
hundreds of thousands of self-gravitating (super-)particles in a stratified equilibrium 
state. The particles were then pushed into a self gravitating regime that eventually led 
to gravitational collapse. Simulations with and without physical collisions were studied 
for a range of particle numbers to test convergence. 

In cases without physical collisions, convergence can not be achieved. However, it 
was possible to change multiple free simulation parameters at the same time, namely 
the particle number N and the smoothing length b, such that in special cases the results 
did not depend on the particle number. We note however that there is a free parameter 
in the simulation (i.e., the smoothing length) that effects the outcome and makes it 
impossible to find the real physical solution. 

In the proto-planetary disc, physical collisions dominate over gravitational scatter- 
ing. In the case of planetesimal formation, the system is marginally collisional. Physical 
collisions will partially randomise the particle distribution, whereas a smoothing length 
does not because the softening length is very large compared to the gravitational cross- 
section. Even if the gravitational particle-particle scattering becomes important, it can 
be seen from equation 15.21 that the velocity dependence of the cross-section ctg differs 
fundamentally from the velocity independent cross-section ac- We therefore argue that 
the super-particle approach should not be used when collisions become important. 

With collisions, the simulation outcome is both quantitatively and qualitatively very 
different from simulations with no physical collisions. The initial size of the clumps is 
larger and the number of clumps is smaller. As there is no smoothing length, there is 
effectively one fewer free parameter. Changing the particle number while keeping the 
collision time constant does not change the outcome. We therefore call these simulations 
converged. 

Additional tests including inelastic collisions with a normal coefficient of restitution 
of 0.25 have been performed to confirm that the results do not depend on the way 
physical collisions are modelled. As expected, the qualitative outcome in terms of 
number of clumps and clump size is different because the physical properties of the 
system have been changed. However, once individual collisions are resolved the results 
converge in exactly the same manner as in those simulations presented above which 
have a coefficient of restitution of 1. 

This chapter focused on the numerical requirements to study gravitational instability 
and the formation of planetesimals in proto-planetary discs. The initial conditions were 
chosen such that the equilibrium is a well defined starting point for the convergence 
study. The collision, damping, collapse and orbital timescales are all of the same order, 
which is effectively the worst case scenario. To provide any constraints on planetesimal 
formation itself, one would have to properly simulate the turbulent background gas 
dynamics which is a far bigger task. 
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However, we were able to determine the numerical resolution needed to resolve dust 
particles in a proto-planetary disc. Assuming that one wishes to simulate dust particles 
realistically, including collisions, and that the particles are uniformly distributed in a 
box of base and height H, the size of each particle is determined by the collision 
rate in the real disc (see equation I5.4p 



where a and are the size and number of particles in the volume L^H of the real disc. 
The number of particles in the simulation A'num has to be large enough to ensure that 
the filling factor is low until clumps form. The requirement that the filling factor is less 
than unity at t = is given by 



where is a safety factor. Although C/ < 1 would be enough to resolve collisions 
initially, it is insufficient to simulate the collapse phase. During the collapse, the filling 
factor rises rapidly. Assuming that one wishes to resolve collisions correctly when the 
particles have contracted by one order of magnitude, one has to include a safety factor 
of Cf = 10~^. If the collapse occurs mainly in two dimensions, as in the simulations 
presented here, a factor of Cf = 10^^ is sufficient. Furthermore, one has to ensure that 
the box contains at least a few unstable modes (see equation I4.29p . All this together 
places tight numerical constraints on numerical simulations of planetesimal formation 
via gravitational instability. 

In the future, the discussion should be expanded to include more physics, namely 
a proper treatment and feedback of the background gas turbulence. This will allow us 
to further clarify the behaviour of planetesimals in the early stages of planet formation 
and eventually test the different formation scenarios. 




(5.7) 




(5.8) 
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Chapter 



Summary 



If I have seen further it is only by standing on the 
shoulders of giants. 

Issac Newton, in a letter to Robert Hooke, 1676 

Several hundred extra-solar planets have been discovered in recent years. This led to a 
revival of planet formation theory. The remarkable differences of extra-solar planetary 
systems and our own solar system keep challenging long standing theories. 

In this thesis, four main results are presented and summarised below. They can be 
categorised into two groups. One group is concerned with the formation of resonant 
multi-planetary systems. The other group discusses the effects of stochastic forcing on 
orbiting bodies such as planets, moonlets and dust. 

6.1 Main results 

6.1.1 Most resonant planetary systems are stable 

In chapter 131 we develop a consistent framework for the treatment of stochastic forces in 
celestial mechanics. The growth rates of all orbital parameters are given as a function 
of the diffusion coefficient in the most generic way. Furthermore, several important cor- 
rection factors which take into account a finite correlation time are included. Previous 
studies have overestimated the diffusion efficiency by more than an order of magnitude 
by ignoring these effects. 

We apply this framework to a single, stochastically forced planet and to two planets 
that are in resonance. The results show that most planetary systems are stable against 
the random forces originating from MRI turbulence. This remains true even when we 
take into account the considerable uncertainties in the amplitude of those forces. 

We find that resonant planetary systems have two different modes, which respond 
differently to stochastic forcing. Thus, systems that did undergo a random walk phase 
during their evolution should show evidence of that. One system, in which this might 
have been observed, is HD128311. We present several successful formation scenarios 
of this system which can reproduce the observed libration pattern. This cannot be 
obtained using laminar convergent migration alone. 



107 



The framework has a wide range of apphcations including Saturn's rings (see chapter 
H]) and promises to have a large impact on future studies involving stochastic forces of 
any type. 



6.1.2 First prediction of orbital parameters of exo-planets 

The planetary system HD45364 is in a 3:2 resonance. The system cannot be formed by 
convergent migration with moderate migration rates. The final configuration is always 
a 2:1 resonance, not the observed 3:2 resonance. 

In chapter [21 we present a successful formation scenario that involves a short phase 
of rapid inwards type III migration. A disc with a mass about five times that of the 
minimum mass solar nebula is sufficient to drive type III migration. Furthermore, a 
large disc scale height is favourable for the long term sustainability of the resonance. 
These results imply that planets form in discs which are more massive than previously 
thought. 

Finally, the outcome of hydrodynamical simulations are directly compared to the 
measured radial velocity curve. We find a new orbital solution, based on the simulation 
results, that has a lower value than the previously reported best fit. The eccentricities 
of both planets are considerably smaller. The results are generic, as in all simulations 
we perform and that form a 3:2 commensurability, the orbital parameters are almost 
identical. This difference is the first prediction of the precise orbital configuration of an 
extra-solar planetary system. It can be verified with new radial velocity measurements 
within the next few years. 

This study opens a new window to quantitative and predictive studies of the dy- 
namics of planetary systems beyond our own solar system. 



6.1.3 Moonlets in Saturn's rings migrate stochastically 

An analytic estimate of both the random walk in semi-major axis and the mean equi- 
librium eccentricity of a small moonlet which is embedded in Saturn's rings is presented 
in chapter |H We show that the most important effects that determine the dynamical 
evolution of the moonlet are collisions with ring particles. If the surface density of the 
rings is sufficiently large that gravitational wakes can form, then those exert an addi- 
tional force on the moonlet which becomes dominant eventually. No matter how the 
moonlet is excited, it will migrate stochastically, with its semi-major axis changing on 
average like y/t for large t. 

We perform realistic three dimensional collisional N-body simulation of up to half a 
million self-gravitating ring particles and confirm the analytic estimates within a factor 
of two. New, pseudo shear periodic, boundary conditions are used which allow us to 
use a much smaller shearing patch than previous authors and we therefore get a large 
performance boost. 

The random walk in semi-major axis leads to changes in the longitude that are 
observable with the Cassini spacecraft. If confirmed, this can be used to calculate an 
upper limit on the age of the propellers and eventually the ring system itself. 
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6.1.4 Physical collisions are vital in numerical simulations of 
planet esimal formation 

In chapter we revisit the numerical requirements in simulations of planetesimal forma- 
tion via gravitational instability. We find that the inclusion of collisions is indispensable 
to obtain a consistent model. If collisions are not included, gravitational scattering is 
the only small scale process present. It will therefore become very important as soon as 
clumping occurs. Naively, one might think that gravitational close encounters could be 
treated as if these were physical collisions. However, properties such as the cross section 
are different compared to physical collisions and gravitational scattering can therefore 
not resemble the real system. 

The results show that earlier studies should not be trusted beyond the initial clump- 
ing phase. Estimates on the minimum number of particles required for convergent 
results are given. 

Future studies on planetesimal formation will have to address this issue and the 
method presented here shows a clear way on how to do that. This will lead to self- 
consistent results in future simulations. 



6.2 Future work 

6.2.1 Multi-planetary systems 

The results on multi-planetary systems in this thesis, as successful as they might be, 
have been restricted to only two systems. However, over a third of the discovered 
planets are confirmed to be in multi-planetary systems. 

Having developed all necessary tools in chapters [2] and [3l it is now an easy task to 
study many systems in a similar fashion. In the following years, I will select up to 15 
multi-planetary systems and study their evolution with both N-body as well as two and 
three dimensional hydrodynamic simulations. This will be the first survey of realistic 
formation scenarios among multi-planetary systems. The results will be more precise 
and not as degenerate as population synthesis models of single planets. 

New observations suggest that a substantial fraction of hot Jupiters might be in ret- 
rograde orbits (Triaud et al. 2010, in preparation). There are vari ous ways to form these 



syste r as, such as Kozai Mig ration or planet-planet scattering fiFabrvcky &: Tremaine 



20071 : iNagasawa et al.l 120081 ) . So far, these methods have not been studied in combi- 
nation with the traditional ideas of planet-disc interaction. I will use Prometheus to 
study highly inclined systems in the presence of gas discs. 

I also want to continue studying the effects of MRI turbulence on planetary systems. 
The strength of stochastic forces coming from turbulence is not well constrained yet. 
For that reason a large interval of possible values of the diffusion coefficient had to be 
studied in chapter [3J Using the MHD module of Prometheus or a different MHD code 
such as Athena, I will explore the parameter space in a systematic way with global non- 
ideal MHD simulations. The measurements of statistical properties of MRI turbulence 
will improve the predictive capability of our analytic random walk model significantly. 

One needs at least 16 (maybe 32 or even 64) grid cells per scale-height to simulate 
the turbulent proto-stellar disc accurately. Assuming one wants to simulate a radial 
domain ranging from 1 to 4 AU and 6 scale-heights in the vertical direction, one needs a 
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96x480x502 grid. This is a challenging task. During my PhD, I have gained knowledge 
on how to speed up calculations using CPU specific hardware features. If the perfor- 
mance of the code becomes an issue and it turns out that higher resolution is required 
to simulate the relevant scales in the system, I can implement these changes to speed 
up the simulations. 

With these new results, I am going to show that a broad selection of systems can 
be formed in a turbulent disc and prove if there is statistically significant observational 
evidence of turbulence in the exo-planet distribution or not. 



6.2.2 Saturn's rings 

The work presented in chapter H] was limited to a rather small subset of the parameter 
space. I plan to extend this study to a wider range of parameters, especially to higher 
internal densities and a broader size distribution to reduce granularity. As soon as 
observational data has been published, I will investigate the constraints arising from 
the measured random walk of moonlets. If the random walk is indeed as strong as 
predicted, it is possible to obtain an upper limit on the age of the ring system. 

Direct simulations of planetary rings are constrained to a local model. Although 
this is a valid approximation in most circumstances, it has also certain drawbacks. For 
example, simulations presented in chapter H] capture the rand om walk in semi-m a i or 



axis but do not capture the net migration over long timescales fICrida et al.ll2010f ). 

I plan to extend direct simulations of Saturn's rings to a global model. Periodic 
boundary conditions in the azimuthal and pseudo shear periodic boundary conditions 
in the radial direction, in a similar manner to those presented in chapter [Hand appendix 
[Pt can be used to extend the shearing box to a wedge that extends radially over several 
tens or hundres of kilometres. We can get rid of Hill's equations and solve the planet's 
gravitational force directly. 

Numerically, this is a challenging task, as the number of particles will be very large. 
Therefore, the particles have to be distributed across processors. As outlined in some 
detail in appendix [Dl MPI is not ideally suited for parallelising tree structures. In this 
case, however, the fact that one dimension is much larger than the others comes in 
handy because boxes far away do not have to communicate directly with each other so 
that their interaction can be pre-calculated. 



6.2.3 Planet esimal formation 

The results on planetesimal formation presented in this thesis are preliminary, showing 
that it is indeed possible to simulate the gravitational collapse of a self-gravitating 
planetesimal population. To do this self consistently, one has to include collisions and 
gas dynamics. The former part has been completed in chapter O The latter part has 



been done by other authors (e.g. iJohansen et al.l 120071 ). A combination of both is the 
logical next step. 

Many groups are currently working on planetesimal formation in a turbulent disc. I 
plan to focus on the differences between global simulations and the local shearing box 
approximation which has been used by most authors. This is essential to understand the 
importance of effects such as the streaming instability which rely on a large over-density 
in the solid to gas ratio. In a local simulation with periodic boundary conditions, this 
ratio is set by the initial conditions and cannot change. However, in a real proto-stellar 
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disc, as well as in global models, the ratio is controlled by large scale structures in the 
disc such as pressure and temperature gradients. 

To gain a better understanding of these global effects, I will couple the new method 
that treats particles correctly to a global MHD simulation. The method to solve the 
non-ideal MHD equations will be the same as for the study of multi-planetary systems. 
To simulate the self- gravity and collisions in the dust layer, I will re-use the GravTree 
code, presented in appendix |Dl which is ideally suited for this project and already fully 
developed and tested. 

6.2.4 Numerical codes 

The primary codes that have been written for this thesis, Prometheus and GravTree, will 

be made publicly available. Both codes will be hosted on github at http : // g ithub . com/hannorein/ 

under the open source license GPL. They are easily portable to different architectures 

as they are written in C and C++, respectively, using only standard hbraries such as 

OpenGL, OpenMP and libpng. Compilation is controlled by the GNU autotool suite. 

Providing the source code is, in my opinion, a necessity as numerical codes are part 
of the initial conditions. All results should be easily reproducible. Firstly, this helps to 
make progress in the subject, and secondly, it also increases the chance to find bugs in 
the code. It is therefore a win-win situation. 



Ill 



Appendices 
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Appendix 



Orbital elements 



In this appendix, the orbital elements of an elliptical orbit are defined. It is a solution 
to the two body problem, the simplest non-trivial problem in celestial mechanics. 

We consider a test particle on an orbit around a central object with mass m at the 
origin. The specific force on the particle is given by F = — ^r, where r is the position 
of the particle. Thus, the equation of motion is simply 

-r. (A.l) 



Taking the vector product of the above equation with r evaluates to zero on the right 
hand side and gives after integrating 

r X r = h, (A.2) 

where h is a constant of motion. The absolute value of h is the specific angular mo- 
mentum. The position and velocity vectors are in the same plane and allowing us to 
restrict ourselves to the two dimensional case. By switching to polar coordinates r and 
(p, we can rewrite equation \A.2\ giving h = r'^cj). We also rewrite equation IA.lt giving 



Gm 



- (A.3) 



The solution to equation IA.3I can then be written as 

V 



1 + e cos(0 — w) ' 



(A.4) 



where e is the eccentricity, w is the longitude of pericentre and p the semilatus rectum 
which can be written as p = h'^{G'm)~^. We restrict the further discussion to an elliptic 
orbit for which < e < 1. In that case p = a{l — e^), where a is the semi major axis. 
We define the true anomaly as f = (p — w. We can further more define an average 
angular velocity, n = 27r/T, called the mean motion, where T is the orbital period, 
given by 

= (A.5) 
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For a circular orbit, the mean motion is simply Vc/a, where Vc = \jGm/a. Now that we 
have a mean angular velocity, we can also define a mean angle M (the mean longitude), 
measured from the pericentre of the planet, via 

M = n{t-T). (A.6) 

Here, r is the time of pericentre passage. M grows linearly in time. 

Finally, we would like to get an equation for r(t) and (p(t). The time dependence 
of equation IA.4I is hidden in 0. It turns out that an explicit form does in general 
not exist and the solution can only be calculated iteratively. We can see this by in- 
troducing another angle, the eccentric anomaly E (see figure ETT]) . It can be shown 



flMurray fc Dermottll2000l ) that a relation between the mean longitude and the eccentric 



anomaly exists such that 

M = E - esinE. (A.7) 

This equation is transcendental but can be solved iteratively, for example with a 
Newton-Raphson method. 

Most of the quantities defined above are plotted in figure lA.ll for illustrative pur- 
poses. 
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Figure A.l: Orbital elements of an elliptical orbit. 



Appendix 



Dissipative forces 



Additional terms can be added to the equations of motion in N-body simulations to 
describe eccentricity and semi-major axis damping. These forces could result from 
various sources, incl u ding a gaseous disc, a planetesimal disc or ti des. The derivation 
follow s iLee fc Peald ( 1200 ll ) and some expressions are taken from iMurrav fc Dermott 
(l2000h . 

We have two additional, non-conservative terms in the equations for both the velocity 
and acceleration. In a Cartesian coordinate system, these are 



dx 
It 

dx 
Itt 



+ 



+ 



dx 

dx 
'di 



dx . dx . 

da de 



dx . dx . 

da de 



and 



(B.l) 
(B.2) 



The expressions are analogue for the y and z coordinates. In order to differentiate 
with respect to a and e, we have to write the cartesian coordinates in terms of orbital 
elements. The positions can be expressed as 



X = r cos cos (a; + /) — r cosz sin sin (w + /), 
y = r sin cos {u + f) + r cos i cos fl sin {u + /) 
z = r sini sin (w + /), 



and 



(B.3) 



which simply follow from the definition of the orbital parameters (see appendix R|) . By 
differentiating once, we get expressions for the velocities 



X 



cosQ[r cos (w + /) — rf sin (w + /)] 

— cosi sinfi[f sin (w + /) + rf cos {u + /)], 

sin fi[r cos (cj + /) — rf sin {u + /)] 

-fcoszcosfi[r sin (w + /) + r/cos {u + /)], 

sin i[r sin {u + f) + rf cos {u + /)]. 



(B.4) 



We also need expressions for the radius r, the time derivative of the radius f, and rf 



119 



in terms of a, e, and /: 



r ^ ^ (B.5, 

1 + e cos / 

rfe sin f 

r = — B.6 

1 + e cos / 
Tin 

rf = ^==(l + ecos/). (B.7) 



e 



To complete the list of equations, we compute several other derivatives needed to cal- 
culate the partial derivatives in equation IB.ll and IB.2I 



dr 
da 

dr 
de 

dr 
da 

dr 
de 

djrf) 
da 

djrf) 
de 



r 
a 



2er cos / 



1 - e2 a(l - e2 



r 

2^' 



r 



e(l-e2) 



2a' 

r/(e + cos/) 
(1 -e2)(l + ecos/)' 



(B.8) 



Putting everything together, we arrive at the following expressions 



dx 


dx 

a 


Itt 


dy 


dy 


dt 


dt 

a 


dz 


dz 


'dt 



x 

—a + 
a 



—a + 
a 



z 

—a + 
a 



r 




l + e2 


a{l - 


e2) 


1 -e2 


r 




l + e2 


a{l - 


e2) 


l-e2 


r 




l + e2 


a{l - 


e2) 


l-e2 



-e, 
e 



?/ . 

-e, 

e 



z . 
-e. 

e 



(B.9) 
(B.IO) 
(B.ll) 
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The additional terms for each of dx/dt, dy/dt, dz/dt are distinct for variations in e: 

dx 



dx 
It 



+ 



dt 



dy 
dt 



dt 



dz 



+ 



dz 
'dt 



X 



a + cos Vl 

2a 



cos i sin Q 



— a + smil 
2a 



+ cos i cos Q 



df d(rf) 
— cos (u + f)- sin {u + f) 



f sin(c. + /) + ^cos(a; + /) 



e (B.12) 



dr 

— cos {u + /) 



d{rf) 
de 



sin {u + f) 



e (B.13) 



|sin(c. + /) + ^cos(a; + /) 



e, 



2a 



d + sin i 



|sin(c. + /) + ^cos(c. + /) 



e. (B.14) 



The final equations IB.9tlB.14l are then used to calculate the dissipative forces and in- 
cluded in the N-body code (see section [2^ . Note that the orbital elements have to be 
calculated every sub-timestep. Depending on the integrator and the required number 
of force evaluations, this can be numerically expensive. 
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Appendix 



Response calculation of particles on 

horseshoe orbits 



In this appendix, we calculate the amplitude of epicyclic motion that is induced by an 
eccentric moonlet in a small particle on a horseshoe orbit. 



C.l Interaction potential 

Due to some finite eccentricity, the moonlet undergoes a small oscillation about the 
origin in the local Hill Cartesian coordinate system (see section 14.1. ip . Its Cartesian 
coordinates then become {X, Y), with these being considered small in magnitude. The 
components of the equation of motion for a ring particle with Cartesian coordinates 
[x, y) = ri are 

— -2fi- = Sr^x-— ^ and (C.l) 

^1 + 24 = -i^. (C.2) 

dt^ at mi oy 

where the interaction gravitational potential due to the moonlet is 

^ ^ Gmim2 

^•^ (r2 + i?2_2ri?cos(0-$))i/2 ^ ' ^ 

with the cylindrical coordinates of the particle and moonlet being (r, 0) and {R, $) 
respectively. This may be expanded correct to first order in R/r in the form 

Gmim2 G'mim2-Rcos(0 — $) m a\ 

1,2 = • (^-4) 

The moonlet undergoes small amplitude epicyclic oscillations such that X = £2 cos (fit + 
e),y = — 26^2sin(f2t + e) where e is its small eccentricity and e is an arbitrary phase. 
Then '^1^2 Kiay be written as 

^^^^__Gm^_^,^_ (C.5) 
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where 



, Gmim2r S2{x cos{Qt + e) - 2y sm{Qt + e)) ^. 
*i,2 = -3 (C-6) 

gives the part of the lowest order interaction potential associated with the eccentricity 
of the moonlet. 

We here view the interaction of a ring particle with the moonlet as involving two 
components. The first, due to the first term on the right hand side of equation IC.4I 
operates when 82 = and results in standard horseshoe orbits for ring particles induced 
by a moonlet in circular orbit. The second term, equation IC.6[ perturbs this motion 
when S2 is small. We now consider the response of a ring particle undergoing horseshoe 
motion to this perturbation. In doing so we make the approximation that the variation 
of the leading order potential equation IC. 41 due to the induced ring particle perturbations 
may be neglected. This is justified by the fact that the response induces epicyclic 
oscillations of the particle which are governed by the dominant central potential. 



C.2 Response calculation 

Setting X ^ X + C,x,y y + C,y, where ^ is the small response displacement induced by 
equation IC.6I and linearising equations \C.2\ we obtain the following equations for the 
components of ^ 

§i-2a'Jf . - and (C.7) 

at nil ox 

'§,2n'-k. (C.8) 

at'' at mi oy 

From these we find a single equation for in the form 



dt"^ mi dx J mi dy 

When performing the time integral on the right hand side of the above, as we are 
concerned with a potentially resonant epicyclic response, we retain only the oscillating 
part. 

The solution to equation IC.9I which is such that $, vanishes in the distant past 
[t = —00) when the particle is far from the moonlet may be written as 

Ccc = a cos(fit) + (3 sm{nt), (C.IO) 

where 

/•* Fsm{nt) ^ ^ ^ f' Fcos(fit) , 

a = - -y — -dt and /3 = / — -dt. (C.ll) 

After the particle has had its closest approach to the moonlet and moves to a large 
distance from it it will have an epicyclic oscillation with amplitude and phase determined 
by 

[°° Fsm(nt) ^ , ^ P Fcos(fit) , 
«oo = - / ^] /3oo = / qT^'^^- 
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In evaluating the above we note that, ahhough F vanishes when the particle is distant 
from the moonlet at large it also oscillates with the epi cyclic angular frequency Q 
which results in a definite non zero contribution. This is the action of the co-orbital 
resonance. It is amplified by the fact that the encounter of the particle with the moonlet 
in general occurs on a horseshoe libration time scale which is much longer than Q~^. 
We now consider this unperturbed motion of the ring particles. 



C.3 Unperturbed horseshoe motion 

The equations governing the unperturbed horseshoe motion are equations IC.2I with 

4,,, = *;^ = -^!!!l!!!2. (c,i3) 



We assume that this motion is such that x varies on a time scale much longer than 
so that we may approximate the first of these equations as a; = — [2/(3f2)] {dy/dt). 
Consistent with this we also neglect x in comparison to y in equation IC.13I From the 
second equation we than find that 



d^y 3 d^l 



dt"^ mi dy 
which has a first integral that may be written 



2, (C.14) 



dt J \y\ 

where as before h is the impact parameter, or the constant value of |a;| at large distances 
from the moonlet. The value of y for which the horseshoe turns is then given by 
y = Vo = 2AGm2/{9Q'^b'^). At this point we comment that we are free to choose the 
origin of time t = such that it coincides with the closest approach at which y = y^. 
Then in the approximation we have adopted, the horseshoe motion is such that ?/ is a 
symmetric function of t. 



C.4 Evaluation of the induced epicychc amphtude 



Using equations IC.6| IC.9I and IC.15I we may evaluate the epicyclic amphtudes from 
equation IC.12I In particular we find 



Gm2 £2 coiiiytt 



5 



(3x2 ^ i2y'^ 



(C.16) 



When evaluating equation IC.12| consistent with our assumption that the epicyclic os- 
cillations are fast, we average over an orbital period assuming that other quantities in 
the integrands are fixed, and use equation IC. 151 to express the integrals with respect to 



t as integrals with respect to y. We then find 
where 



— A'ooSine and /3oo = — v4'ooCose 



A. 



y/QGm^£2 



yo 



6fir3^1 - yo/y 



8Gm2 fl 1 



- - - ) + 12|/2 

y Vo 



,-2 



dy, (C.17) 
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with = 8(^7712 /(3fi^?/o)(l — Uo/y) + u"^- From this we may write 



where £if is the final epicychc motion of the ring particle and the dimensionless integral 
X is given by 

We remark that the dimensionless quantity rj = 8G'm2/(3f2^?/o) = 8(r/^'/?/o)^, with 
rn being the Hill radius of the moonlet. It is related to the impact parameter b by 
7] = 2~^{b/rHY ■ Thus for an impact parameter amounting to a few Hill radii, in a very 
approximate sense, rj is of order unity, yo is of order th and the induced eccentricity 
£if is of order £2. 



C.5 Numerical verification 

We verify the analytic calculation with the numerical code GravTree (see appendix iDl). 
A moonlet with Hill radius th = ?)2m. is setup on an eccentric orbit with £2 = 130m 
and Q = 0.0001314 s~^. A ring particle is placed on an initially circular orbit far away 
from the moonlet with different impact parameters b that all lead to horseshoe orbits. 

The amplitude of epicyclic motion of the ring particle, £1, is plotted in figure [QTl 
for four test cases with b = 0.25rH, b = O.drn, b = O.Tbrn and b = th- We also plot 
the analytic estimates of the final epicyclic motion given by equation IC.18I as straight 
lines. The agreement is very good, showing some differences at large b. Note that an 
impact parameter greater than Ith ~ 2rii leads to impacts so that the departure from 
the analytic estimate does not influence the circularisation time calculated in chapter 
m significantly. 
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Figure C.l: Epicyclic amplitude Si of a ring particle on a horseshoe orbit as a function 
of time for different impact parameters b. The horizontal lines are the estimate given 
by equation IC.18I 
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Appendix 



N-body simulations with GravTree 



This appendix describes the GravTree code that has been written to solve many body 
systems. The particles might interact gravitationally and through physical collisions. 
Modules incorporating gas drag and stochastic forces have also been implemented. 

First, the relevant equations of motions are presented, followed by a description of 
the special symplectic integrator used for these. Second, the Barnes Hut tree code is 
explained. Finally, tests and optimisation strategies that have been implemented are 
discussed. 



D.l Hill's equations 

In many situations it is impossible to simulate an entire accretion disc or planetary 
ring. For example, in the case of a proto-planetary disc there are about one trillion 
(10^^) metre-sized particles. Numerically, we can only simulate a few million, at most. 

Fortunately, it is usually sufficient to simulate a small patch of the disc that in- 
cludes most of the physics. In such a patch (also known as a shearing sheet) the shear 
is linearised and geometric effects are ignored. This leads to a well defined set of equa- 
tions with shear periodic b oundary conditions in the radial and azimuthal directions 



(IWisdom &: Tremainelll988l ) 



Let us start from the equations of motion of a test particle around a single star at 
the origin with mass M such that the gravitational force can be described by 

GM , , 

r = (D.l) 

where r is the distance from the star. We now move to a non-inertial frame, rotating 
with angular frequency f2 = |r2|. The time derivative of a vector quantity v in the new 
frame is given by 

Vr = Vi - ri X V, (D.2) 

where subscripts i and r denote the inertial and rotating frames, respectively. The 
acceleration in the rotating frame is thus given by 

ar = ai — 2 f2 X Vr — ri X (fi X Tr), (D.3) 
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where is the velocity and is the position vector in the rotating frame. Adopting a 
coordinate system in the rotating frame where x corresponds to the radial direction, y 
to the azimuthal direction and z to the rotation axis, setting the origin at a distance a 
away from the star, we can simplify the above equation to get 

X = ri + 2n y, + n"^ {x + a)r and (D.4) 
y = -2n Xr + Vr- (D.5) 

Linearising equation ID.ll around r = a + x gives 

GM GM GM GM , , 

+ 2— -a; + ... (D.6) 



(a + xY 

= -n'^a + 2n^x + ... (D.7) 

where we have assumed a Keplerian rotation, f2^r^ = GM. Thus, the linearised equa- 
tions of motion for an unperturbed particle orbiting in the rotating frame are 

X = 3Q^x + 2Qy and (D.8) 
y = -2nx, (D.9) 

where we have dropped the subscript r. We can also linearise the gravity in the vertical 
direction to get 

z = -n'^ z. (D.IO) 

This set of local approximations is known as Hill's equations. We solve these with a 
symplectic leap frog, kick drift kick, time-stepping scheme which is discussed in more 
detail in the next section. 



D.2 Symplectic integrator 

Symplectic integrators in celestial mechanics have several attractive features, such as 
conservation of energy up to machine precision. We would therefore like to use such an 
integrator for solving Hill's equations numerically. A widely used standard symplectic 
integrator is leap-frog which can be described by a kick, followed by a drift and another 
kick sub-time-step (KDK). A possible C implementation of one full time-step with 
length dt is 



void kdk (double* 




double* V, double* a, double 


// Kick 






v[0] += 0.5 * 


dt 


* a[0] ; 


v[l] += 0.5 * 


dt 


* a[l] ; 


v[2] += 0.5 * 


dt 


* a[2] ; 



// Drift 

x[0] += dt * v[0] 

x[l] += dt * v[l] 

x[2] += dt * v[2] 
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// Kick 

update_f orce (x , v , a) ; 
v[0] += 0.5 * dt * a[0] 
v[l] += 0.5 * dt * a[l] 
v[2] += 0.5 * dt * a [2] 



Here, x, v and a are pointers to the position, velocity and acceleration vectors, respec- 
tively. This integrator is symplectic if the forces calculated in update_f orce(x,v,a) do 
not depend on the velocities. However, Hill's equations do depend on the velocities (see 
equations ID.8I and ID.9|) . Note that the z component (equation ID.10|) doe s not depend 



on th e velocities and thus can be integrated with a standard method. iQuinn et al 
( I2OIOI ) describe a new symplectic integrator for Hill's equations. The changes to the 
KDK algorithm are minimal and only one additional value per particle has to be saved. 
A modified version of the above code is 



void kdk_quinn (double* x, double* v, double* a, double* Py, double dt){ 
// Kick 

v[0] += -0.5*dt * (OMEGA*OMEGA * r [0] - a[0]); 

Py = v[l] + 2.0*0MEGA * r[0] + 0.5*dt * a[l] ; 

v[0] += dt * OMEGA * Py; 

v[l] = Py - OMEGA * r[0] - OMEGA * (r[0] + dt * v[0]); 

v[2] += 0.5*dt * a[2] ; 

// Drift 

X [0] += dt * V [0] ; 

x[l] += dt * v[l] ; 

X [2] += dt * V [2] ; 

// Kick 

update_f orce (x , a) ; 

v[0] += dt * OMEGA * Py; 

v[0] += -0.5*dt * (OMEGA*OMEGA * r [0] - a[0]); 
v[l] = Py - 2.0*0MEGA * r[0] + 0.5*dt * a[l] ; 
v[2] += 0.5*dt * a[2] ; 



Here, OMEGA is the orbital frequency and Py is initialised to v [1] + 2 . 0*OMEGA * r [0] . 
Py is a conserved quantity for an unperturbed orbit. Note that the drift step is unmod- 
ified. This allows for an unchanged collision detection algorithm, usually performed 
during the drift sub-time-step. The function update_f orceO does not include the 
terms due to Hill's equation anymore, only additional forces (e.g. gravity from other 
particles, gas drag, . . . ). 

The symplectic integrator has proven to be very useful for integrations of moonlets 
in Saturn's rings, where the eccentricity is generally very small (~ 10"®) and the system 
is integrated for hundreds of orbits. 
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D.3 Gravity calculation 



In the iV-body problem that we are considering here, we have to solve Newton's equa- 
tions of universal gravitation for a large number of particles N. The gravitational force 
on the i-th particle is given by 



F.r.. = j:G-^^e.,, (D.ll) 

where h is the smoothing length used to avoid divergencies in numerical simulations. In 
simulations where particles have a finite size and physical collisions are included we can 
safely set 6 = 0. The unit vector in the direction of the gravitational force between the 
i-th and j-th particle is Bij. Calculating the gravity for each particle from each other 
particle results in O (N^) operations. To reduce the number of operations, one can use 
different approximat ions to equation ID. Ill Here, we use a Barnes-Hut (BH) tree code 



(iBarnes Sz Hutlll986l ). The BH tree reduces the number of calculations to O (NlogN). 

A tree code in three dimensions works by sub-dividing the original box into eight 
smaller cells with half the length of the original box. This can be done iteratively, by 
adding one particle to the box at a time. The new particle is given to the root box 
first. If the root box has already been sub-divided into smaller boxes, the particle is 
passed down to the next level. If the particle ends up in a box that has not yet been 
sub-divided, but already contains a particle, then this box is divided into smaller boxes 
and the particle that was in this cell is removed. The original and the new particles 
are then re-added. This process is repeated until all particles are inserted. The depth 
of the tree in a homogeneous medium is approximately logg A^. The next step is to 
pre-calculate the total mass and the centre of mass for each cell at each level of the 
tree. 

To get the force acting on a particle, one starts at the top of the tree and descends 
into the tree as far as necessary to achieve a given accuracy. If the current cell is far away 
from the particle for which the force is calculated, then the detailed density structure 
within this cell is not important. All that matters is the box's monopole moment (total 
mass and centre of mass)^. One therefore does not have to descend into this branch 
any further. This criteria can be quantified by an opening-angle 6, which is defined as 
the ratio of the box width and the distance from the particle to the centre of the box, 
as shown in figure [PTTl In all simulations we use 6 < 0.5. 



D.4 Boundary conditions 

Our implementation can handle arbitrary aspect ratios of the box, which is especially 
useful if a small annulus has to be simulated. We use ghost rings made of ghost boxes 
in the radial and azimuthal directions as shown in figure ID. 21 The figure shows one 
ghost ring. A ghost box is either simply a shifted copy of the main box or a completely 
separate simulation, depending on the type of simulation (see below). 

The gravity on each particle is then calculated by summing over contributions from 
each (ghost) box. This setup approximates a medium of infinite horizontal extent and 
avoids large force discontinuities at the boundaries. In general no ghost boxes are used 
in the vertical direction because the disc is stratified. 

^This could be easily extended to include higher order moments. 
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Figure D.l: The opening angle 6 is the acute angle in the triangle and defined as the 
ratio of cell width divided by the distance from the particle to the cell. If ^ < 0.5, the 
cell is sufficiently large away and any substructure can be ignored. 



D.4.1 Shear periodic boundary conditions 

In a standard simulation, we consider the main box together with eight ghost boxes that 
are all identical copies of the main box shifted according to the mean shear. On account 
of this motion ghost boxes are removed when their centres are shifted in azimuth by 
more than 1.5 times their length from the centre of the main box and then reinserted 
in the same orbit on the opposite side of the main box so that the domain under 
consideration is prevented from shearing out. 

If a particle in the main box crosses one of its boundaries, it is reinserted on the 
opposite side so that the total number of particles is conserved. The gravity acting on 
a particle in the main box, is calculated by summing over the particles in the main box 
and all ghost boxes. 

A finite number of ghost rings can only act as an approximation of an infinite 
medium and the gravitational force will tend to concentrate particles horizontally in the 
centre of the box. In simulations of planetesimals (see chapter E]), we use 8 ghost rings, 
corresponding to 63 ghost boxes. The asymmetry of the gravitational force between 
the centre and the faces of the box is reduced by a factor of 20 compared to using no 
ghost rings at all. Note, however, that a small asymmetry remains in every simulations 
and lea ds to higher a conce ntration of particles in the box centre. In other simulations. 



such as iTanga et al.l (l2004j ). this effect can not be seen because the system is initially 
gravitationally unstable and integrated for only one orbit. Some systems that we are 
interested in are marginal gravitationally unstable and this slight asymmetry might 
become important after many orbits (see chapter [5]). 
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Figure D.2: Shearing box and ghost boxes. The ghost-boxes move in the y direction, 
following the mean shear, leading to shear-periodic boundary conditions. Here, the 
boxes have an aspect ratio of 2, but any other aspect ratio is possible. 



D.4.2 Pseudo shear periodic boundary conditions 

If a strong perturber is present in the simulation, such as a moonlet in chapter HI we use 
a different setup that we call pseudo shear periodic boundary conditions. In that case, 
all ghost, or auxiliary, boxes are identical copies whose centres are shifted according to 
the background shear as in a normal shearing sheet. However, the auxiliary boxes are 
a completely separate simulation, mimicking the background state of the ring system 
without a perturber. 

If a particle in the main box crosses one of its boundaries, it is discarded. If a particle 
in the auxiliary box crosses one of its boundaries, it is reinserted on the other side of 
this box, according to normal shear periodic boundary conditions. But in addition it is 
also copied into the corresponding location in the main box. 

Other simulations ( Lewis fc Stewart 20091 ) use a very long box (about 10 times 



longer than the typical box size that we use) to ensure that particles are completely 
randomised between encounters with the moonlet. We are in general not interested in 
the long wavelength response that is created by the moonlet. Effects that are most 
important for the moonlet 's dynamical evolution are found to happen within a few Hill 
radii. Using the pseudo shear periodic boundary conditions, we ensure that incoming 
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Figure D.3: Resolving collisions. Top left: An overlap has been detected in the original 
frame. Top right: Reference frame where one particle is at the origin and at rest. 
Bottom left: Collision is resolved using momentum and energy conservation. Bottom 
right: Back in original frame. 



particles are uncorrelated and do not contain prior information about the perturber. 

The gravity acting on a particle in the main box, which also contains the moonlct, 
is calculated by summing over the particles in the main box and all auxiliary boxes. 
The gravity acting on a particle in the auxiliary box is calculated the standard way, by 
using ghost boxes which are identical copies of the auxiliary box. 

This setup speeds up our calculations by more than an order of magnitude. 



D.5 Collision detection 

The collision detection algorithm makes use of the already existing tree structure to 
search for nearest neighbours. Again, this reduces the complexity from 0{N'^) to 
O(A^logA^), and can easily be done in parallel. We search for overlapping particles 
that are approaching each other. The order in which multiple collisions involving the 
same particle are resolved are not important in our case. We can ignore this effect as 
long as the time-step is sufficiently small. 

This approach might lead to another problem. Consider two particles that are con- 
stantly touching; they might sink into each other over long timescales due to numerical 
errors and other dissipative effects. To avoid this fate, we hmit the minimal impact 
velocity V|,min to a tiny fraction of the random velocity dispersion, vi^min has to be small 
enough to not affect the global evolution of the system. 
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Collisions are then resolved using energy and angular momentum conservation and 
a (velocity dependent) coefficient of restitution. This is done by moving to a reference 
frame in which one particle is at rest at the origin and the other particle is on the 
x-axis, as shown in figure ID. 31 

D.6 Optimisations 

Cache 

The tree structure leads to a large number of random memory accesses. This can 
potentially harm the performance of the code if many cache misses occur. When a 
cache miss occurs, the CPU has to request the data from the main memory, which is 
much slower than a request that can be handled by the cache. Using a simple heuristic 
branch prediction, one can speed up the calculation. We do this by pre-caching lower 
level tree branches before working on the current level. In C++, using non-standard 
MMX extensions, a typical tree walking subroutine looks as follows: 

#include <xmmintrin.h> 
void cell : : walk_tree()-[ 

// Try to cache position and tree structures of next lower level 
for (int i =0; i<8; i++) i 
if (d[i] !=NULL)-[ 

_iniii_pref etch ((const char*) &(d[i] ->pos) , _MM_HINT_NTA) ; 
_min_pref etch ((const char*) d[i]->d, _MM_HINT_NTA) ; 

} 

} 

//Do calculations here 
// ... 

// Then descent into next level 
for (int i =0; i<8; i++){ 
if (d[i] !=NULL)-[ 

d [i] ->walk_tree ; 

} 

} 

} 

While calculations on the current level are performed, the CPU can try to cache data 
for future calculations in idle CPU cycles. This trick can give a performance boost of 
up to 20%. 

MPI 

Parallelising tree codes on distributed memory systems is hard. This is not a shortcom- 
ing of the algorithm itself but rather the implementation of interconnect libraries such 
as MPI. A standard MPI program knows in advance when to send and receive data. 
This is a reasonable approach for grid based hydrodynamic codes for example, where 
exactly the same boundary cells have to be transmitted each time-step. 
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In a tree code, as described above, the situation is different. It is not possible to 
predict which cells will have to be opened during the calculation of the gravitational 
force for one particle. Thus, if the cell structures are distributed over many machines, 
we have to pull the required data instantaneously. 

With the current implementations of the MPI library, this is not possible without 
large overhead. Thus, we use a brute force method to do calculations on distributed 
memory machines. 

Each machine holds an exact copy of the entire tree. The tree construction and 
reconstruction is done on a single node and only the computationally expensive calcu- 
lations are done in parallel. These include gravity calculation and collision detection. 

This kind of parallelisation is limited by the number of particles, as the memory 
usage of each node scales with A^logA^. Furthermore, the communication (usually 
all-to-all and one-to-all) scale also as A^logA^. Nevertheless, it has proven to be very 
practical up to ~ 10^ - 10^ 



D.7 Tests 
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Figure D.4: Conserved quantities in the shearing box. The top panel shows the centre 
of epicyclic motion relative to the centre of the box. The bottom panel shows the 
amplitude of epicychc motion. 



We begin with testing several conserved quantities. In figure [Dl4l the amplitude and 
the centre of epicyclic motion of a test particle in the shearing box are plotted as a 
function of time. One can see that both quantities are conserved accurately without 
any asymmetry over many dynamical timescales, as expected for the symplectic inte- 
grator. Furthermore, we made sure that the treatment of collisions is working correctly, 
especially across ghost cell boundaries. 
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Next, we present a comparison between GravTree and the results obtain by lSchmidt et al. 



(120091 ). This is a comprehensive sample of runs including of both self-gravity and col- 
lisions. In every simulation, the particle density is p = 0.9g/cm'^, and the box width 
corresponds to four critical Toomre wavelengths. The coefficient of restitution is 0.5. 
By changing the optical depth and semi-major axis of the planetary ring, one can sweep 
through different regimes. Snapshots of the den sity distribution are plotted in figure 
ID. 51 Those should be compared to figure 14.7 in lSchmidt et al.l ( 120091 ). 

One can clearly see different regimes. On the right hand side of the plot, particles 
clump into aggregates. On the top right boundary, strong gravitational wakes are 
visible. For r = 1.8, a = 70000 km and simulations close by, viscous ov erst ability can 
be obs erved as strong azimuthal structures. Overall, the agreement with lSchmidt et al 



( 120091 ) is very good. For large a, the simulations performed with GravTree have a trend 
to form aggregates slightly earlier. 



As a further, more quantitative test, we reproduce the results from iDaisaka et al. 



( I2OOII ). Therefore, the translational, collisional and gravitational components of the ring 
viscosity are measured, as defined by equation 19 in said paper. The res ults are shown 



i n figu re ID. 61 and should be compare to the left hand side of figure 7 in iDaisaka et al 
( I2OOII ). The viscosity is in units of the particle radius and the orbital frequency, the 
Hill radius (we here use the slightly different convention = {2mp/ {3M))^^^ a) is 
normalised by the particle radius. One can see that for small Hill radii (low density 
particles), the effects of self-gravity are not important. Again, very good agreement 
with the previous study is achieved. 



D.8 Visualisations 

The visualisation of large data sets plays an increasingly important role in the scientific 
process. GravTree comes with a visualisation module, based on the OpenGL API. It can 
be used to visualise simulations in real time while they are running or in a standalone 
playback mode. The real time mode has proven to be especially useful during the 
development process. 

Beside the particle distribution itself, average quantities such as the velocity disper- 
sion and the filling factor can be shown as a transparent overlay. The tree structure 
can also be plotted. All graphic routines are easily expandable to incorporate any diag- 
nostic that the current project requires. Two screenshots are shown in figures [D?7l and 

El 
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Figure D.5: Dependence of self-gravity wakes on the geometrical optical depth r and the semi major axis a (in units of 1000km). 
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Figure D.6: Translational, collisional and gravitational viscosity components as a func- 
tion of normalised Hill radius of the ring particles. 
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Figure D.7: Screenshot of the OpenGL visualisation module of GravTree, showing a 
simulation of a moonlet embedded in Saturn's rings. The entire computational domain 
is shown with an overlay, showing the Toomre Q parameter as coloured squares. 
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Figure D.8: Screenshot of the OpenGL visualisation module of GravTree, showing a 
simulation of a moonlet embedded in Saturn's rings. A closeup of the moonlet with 
particles resting on the moonlet 's surface (blue) is shown. 
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Appendix 



Hydrodynamical simluations with 

Prometheus 



Simulations of planets embedded in an accretion disc have been performed by a large 
number of authors in recent years. There are two main types of algori thm that can be 



used to simulate g as dynamics: finite diff erence method s such as ZEU S (jStone fc Norman 



19921), NIRVANA fIZieder fc Yorkelll997f) and FARG O (lMassetll2000D or Godunov meth - 



ods like RODEO flPaardekooper fc Mellemal l2nn6h and ATHENA f IStone et all 120081 ). 
As the name suggests, finite difference codes evaluate partial derivatives with a finite 
difference approximation. Godunov methods solve (exactly or approximately) a Rie- 
mann problem at every cell interface. 

In this appendix, the code Prometheus is presented, which is a finite difference code 
in spherical coordinates, specialised for planet disc interactions. It is similar to the 
FARGO code but can also run simulations in three dimensions and has a module for 
solving the full MHD equations. This appendix can not be a comprehensive discussion of 
finite difference codes in general and merely lists the main concepts and implementation 
specific features of Prometheus. 

Prometheus is written entirely in C and requires no external libraries. It can be com- 
piled with a visualisation module, in which case it has to be linked to OpenGL, GLUT 
and libpng. Prometheus will be made publicly available at http : / / github . com/hannorein/prometheu 



E.l Navier-Stokes equations 

For many astrophysical processes collisions between gas particles are frequent enough 
such that a continuum description can be used. Here, the proto-planetary disc is treated 
as a gas whose motion is described by the Navier-Stokes equations. Although an MHD 
module has already been written for this code, it is not described here, due to current 
lack of testing and it has also not been used in scientific runs yet. We therefore restrict 
this discussion to the classical Navier-Stokes equations without magnetic fields, which 
we list in the following for both two and three dimensions. 
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Figure E.l: Spherical coordinate system used in Prometheus 



E.1.1 Three dimensions 



In spherical coordinates as defined in figure \E.1\ the inviscid Navier- Stokes equations 
in conservative form are 



dp 
di 

dpUr 

~dr 

dprug 
dt 

dpru^ cos 9 
di 



+ V ■ (pu) 

+ V ■ (pUrU) 

+ V ■ (prugu) 

+ V ■ {pru^p cos 9 u) 







ui + ui 



r 



-p tan 6* ul 



dp 

dp 

~ d9 
dp 



(E.l) 
(E.2) 

+ pae r (E.3) 
+ pa^f) r cos 9. (E.4) 



Here, p is the density, u 



, Uif,, ug) the velocity, p the pressure and a 



the acceleration due to the star, planets or other forces (e.g. from a rotating coordinate 
system). The terms on the left hand side describe advection, the terms on the right 
hand side are source terms. Equations IE.HIE.4I are the continuity equation, the radial 
momentum equation, the meridional momentum equation and the angular momentum 
equation, respectively. The system of equations is closed with an equation of state, 
determining the pressure p. All simulations presented in this thesis use an isothermal 
equation of state p = cj.p, where the sound speed Cg might vary with radius. 



E.1.2 Two dimensions 

One can simplify the equations in the previous section by averaging over the 9 compo- 
nent, or z for cylindrical coordinates, to obtain a two dimensional system of equations. 
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The vertically averaged surface density is given as the integral over the density 



p(r, 0, z)dz. 



(E.5) 



Thus, in two dimensional polar coordinates the Navier-Stokes equations can be written 



as 



'dt 

dHUr 

dt 
dt 



+ V ■ (Eu) 
+ V ■ (Sm,.u) 
+ V ■ (SrM0 u) 



= 



r 



dp 

dr 
dp 



+ Sa0 r. 



(E.6) 

(E.7) 
(E.8) 



Here, u is the velocity two- vector {ur,u^), sometime also written as {u,v). Equation 
IE.6I is the two dimensional continuity equation, equation IE. 71 is the radial momentum 
equation and equation IE.8I is the angular momentum equation. 

In general the two dimensional system is a good approximation in planet disc simu- 
lations. In two dimensions, a smoothed planet potential is used to simulate the vertical 
structure of the circum-planetary disc (see equation ID.11|) . The smoothing length b is 
approximately half the local disc scale height. 



E.1.3 Viscosity 

Strong physical viscosity is often used to model the effective viscosity resulting from 
turbulence (see section ri.3.ip . This can be achieved by an additional source term in 
equations lE.lllE^ and lE.7IIE.8l Here, we restrict ourselves to the two dimensior ial case 
in wh ich the additional terms to equations IE. 71 and IE. 81 are given by (see e.g. iMasset 



20021) 



= — ^ — \ — ^7 i^-yj 

r or r 0(p r 

r or r 0(p r 

The components of the symmetric viscous stress tensor r are given by 

2 

Trr = Drr Sz/ Vu, (E-H) 

3 
2 

r<^0 = 2 Sz/ D,^^ - - Sz/ Vu and (E.12) 

Tr<t, = r^r = Dr4>, (E.13) 



where we have used 



D„ = (E.14) 

C*. = + ^ (E.15) 

r d(f) r 

1 / d{u^/r) IduA , . 
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Figure E.2: Discretisation and staggered grid in two dimensions, showing the location 
of the cell and face centred variables. 



Here, u is the shear viscosity and the bul k viscosity has been assume d to be zero, u can 
be expressed in terms of the a parameter fIShakura fc Syunyaevlll973l ) via a = v j (c^/fi), 
with Cs being the local sound speed and 27r/f2 being the orbital period. 



E.2 Discretisation 



To solve equations IE.1IIE.4I or IE.6IIE.8I numerically, we discretise them on a grid. A 
staggered grid, where the velocities are face-centred and the densities are cell-centred 
has been used, as illustrated in the two dimensional case in figure [El2l The discretisation 
is identical to that used in the ZEUS code. A staggered grid has the advantage that 
spacial derivatives are automatically second order accurate, but does also add some 
complexity, especially in three dimensions. 

In the radial direction, the boundary conditions can be chosen to be either reflecting 
or damping. In the azimuthal direction, the boundary conditions are periodic. It is also 
possible to simulate only a wedge of the disc (without any planets). In three dimensions, 
periodic boundary conditions are used in the d direction. 
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Figure E.3: A full time-step performed by Prometheus. Dashed boxes represent optional 
modules. 



E.2.1 Operator splitting 

Prometheus makes use of the operator splitting technique, where the equations are 
solved in small chunks. To do that, a time-step is divided into several sub-time-steps. 
Figure IE. 31 illustrates the main loop in Prometheus, including the additional function 
calls for the MHD module, which are not discussed any further here. 

First, the source terms (right hand side of equations IE. ltlE.4l and IE. 6HE.8P are added. 
Then, the advection terms (left hand side of equations IE.ltlE.4l and IE.6IIE.8p are solved 
using the half updated values after the source step, ignoring all terms on the right 
hand side. We further use dimensional splitting, solving the advection in the radial, 
azimuthal and meridional directions consecutively. 

This method is only a simplified approximation to the real solution, but has proven 
to be more accurate than a single time-step. It is (approximately) second order accurate 
in time. Furthermore, it has the advantage that new physics can be easily incorporated 
as additional sub-time-steps. In this code we make use of this when an MHD mod- 
ule is added. However, operator splitting can also lead to spurious results in special 
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circumstances, as shown in section IE.3.21 



E.2.2 Constrained transport 

The transport step is solved using the integral representation of equations IE.mE.4l 
and IE.6HE.8I to ensure conservation of mass and (angular-)momentum up to machine 
precision. We therefore need to have the fluxes at cell boundaries. There are different 
ways of estimating the cell faced, time centred flux at each interface, whose discussion 
goes beyond the scope of this section. In general, higher order methods are less diffusive 
but computationally more expensive. Both the first order donor cell method and the 
second order van Leer method have been implemented. 

Additionally, a concept called constrained transport is used, where the momentum 
flux calculations are constrained, such that they are consistent with the mass flux. This 
leads to improved conservation properties. It is physically motivated, as all conserved 
quantities are advected by the fluid. 



E.2.3 Splitting the angular momentum equation 

Massed ( 2000 ) proposes a method, called FARGO (Fast Advection in Rotating Gaseous 



Objects), to speed up calculations in which there is a large mean background flow, as 
it is the case in a sheared accretion disc. The idea is to rewrite the advection terms, 
the second column in equations IE.mE.4l and IE.6tlE.8[ by replacing the angular velocity 
Uff, by the sum + U(j, + u'^. We choose -u^ to be the mean angular velocity of the 
annulus, rounded such that Ucf, 6t is an integer multiple of the grid spacing 6(f). is 
the remainder of the rounding operation and thus less or equal than 0.5 6(f)/ 6t. Both 
and are independent of and defined for an annulus with constant r. The residual 
velocity is given by u'^. Thus, the contribution from the component in a generic 
advection term of the form V ■ (a u) can be written as 

d d d 

— {a{u^ + u^ + u'^)) = u^—{a)+u^—{a) + —{au'^), (E.17) 

(a) (b) (c) 

where a can be any quantity. The term (c) is just like the old advection term, but with 
a different velocity, as if we were in a rotating frame. The terms (a) and (b) can be 
interpreted as a simple shift. Because we've chosen such that u^/6t is an integer 
multiple of the grid spacing 6(f>, the term (a) corresponds to shifting the whole ring 
of cells, which is numerically very inexpensive. Because Im'^I ^ |m<^| and \U(j,\ ^ \u^\ 
the time-step constraint arising from the CFL condition allows us now to use much 
larger time-steps and the terms (b) and (c) can be solved using a standard advection 
algorithm. The splitting of the advection term gives an efficiency boost of a factor ~ 10, 
in a typical simulation. 



E.3 Tests 



E.3.1 Shock-tube problems 



Shoc k-tube problems are a standard test in numerical codes (see e.g. IStone fc Norman 



19921). Even though these test are in general not very thorough and easy to pass, we 



148 



in 
c 

0) 



2e-05 



1.8e-05 



1.6e-05 



1.4e-05 



1.2e-05 



1e-05 



simulation + 
analytic result 



I 

I 

V 



>+ 



1 IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII N 



^ lllllllllllllllllllllllllllll 



0.2 0.4 0.6 0.8 1 

X 

Figure E.4: Density profile for a strong sliock-tube problem at time t = 10 with Mach 
number v/c = 2.9. 



present a one dimensional example in Cartesian geometry, mainly for comparison with 
a low Mach number shock in section IE.3.21 

The Mach number of the shock front is v/cg = 2.9. An equidistantly spaced grid 
with 256 grid points, a CFL condition of 0.5 and a sound speed of Cg = 0.05 is used. 
The density profile after t = 10 together with the analytic solution is plotted in figure 
IE. 41 One can see that the artificial viscosity smears out the discontinuity over several 
grid cells. Furthermore, the rarefaction wave results in some undershooting. These 
effects are expected and present in any finite difference code. 

E.3.2 Dispersion error 

Finite difference codes are known to have a numerical dispersion error because they 
do not, contrary to Godunov methods, solve the equations along wave characteristics. 
The fluid equations are solved using operator splitting. Each part, namely the source 
and transport step, are evaluated successively. Propagating sound waves involve both 
operators. The advection of a passive scalar, contrariwise, requires only one operator. 

In most situations, this does not create any problems, as the error associated with the 
operator splitting is small and mostly affects the highest frequency modes which appear 
on the grid scale. However, if there is a discontinuity in the flow, large frequencies are 
present that can lead to a large dispersion error. For supersonic shocks, the artificial 
viscosity takes care of this, by smearing out discontinuities over a few grid cells, as 
shown in section IE. 3. II For low amplitude, subsonic shocks this might not be the 
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Figure E.5: Density profiles for different grid resolutions in the isothermal (c^ = 0.05) 
shock-tube problem at time t = 4. The initial density perturbation is 1%. 

case. The shock behaves like a travelling wave package, where the package includes a 
discontinuity. 

Unfortunately, we are exactly in this situation for a low mass planet being embedded 
in an accretion disc. The wake of the planet contains a discontinuity and its amplitude 
depends on the planet's mass. For sufficiently small masses, the discontinuity becomes 
subsonic and we expect noticeable numerical dispersion to occur. Since the wake is 
stationary in a frame co-rotating with the planet, the result will also be a stationary 
pattern. 

In the following, we will first illustrate the problem with a shock tube test in a 
Cartesian box. We will then show that the same pattern appears in simulations of 
low mass planets and outline possible solutions such as physical viscosity and Godunov 
methods. 

Low Mach number shock-tube tests 

We set up a simple test problem that shows the problematic behaviour. We work in a 
Cartesian grid with resolutions 256x256, 512x512 and 1024x1024 and use an isothermal 
equation of state with sound speed Cg = 0.05 in dimensionless units. The box has 
periodic boundary conditions and we set up a shock-tube problem by increasing the 
density in the left part of the box by 1%. The analytic solution after t = 4 is plotted as 
a solid line in figure [K5l Due to the low amplitude of the initial perturbation, the waves 
are basically sound packages travelling at the speed of sound, i.e. the rarefraction wave 
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is very steep. Large ringing can be observed in the numerical solutions, which are also 
plotted in figure IE.5I This is because the numerical dispersion error is largest for the 
smallest wavelength. As expected for a Gibbs-like phenomenon, the over and under- 
shoots do not die out with increasing resolution. However, they get more concentrated 
near the discontinuity. 

Note that in all simulations, the Neumann and Richtmyer artificial viscosity is 
switched on. The artificial viscosity does not kill the oversh oots as it is quadratic 



in the velocities. This is reasonable approach in strong shocks (IStone fc Normanlll992l . 
section IE.3.ip . However, the velocity perturbations are very small compared to the 
sound speed in this case. 

Low mass planets 

A similar phenomenon can be observed in simulations of planet disc interactions. All 
integrations were performed with Prometheus. The results have been confirmed with 
a variety of codes (FARGO, NIRVANA and RH2D) and do not depend on the imple- 
mentation. 

The phenomenon is stronger for low mass planets, as the shocks that appear in the 
wake are weaker and eventually become subsonic. Therefore, a planet with a very low 
mass of TTLp = 1.26 ■ 10~^Mq = 0.42Mearth is used in these simulations. The cylindrical 
grid is spaced equidistantly in the radial direction and has a resolution of 128x384. We 
plot the surface density contours of one simulation with and without physical viscosit y 



in figure IE. 61 We over-plot the analytic position of the wake (jOgilvie fc Lubowi 120021 ) . 

One can clearly see that there are over-shootings occurring, similar to those discussed 
above. Including a small amount of physical viscosity smears out the phenomenon. In- 
creasing the resolution also helps to reduce the strength of the over-shootings. Numer- 
ical diffusion becomes more important if more grid cells, and therefore also a smaller 
time-step, are used. 

For reasonable resolutions (A^^ > 512), the torque on the planet is not significantly 
affected by this effect for two reasons. Firstly, the over-shootings appear a few scale 
heights away from the planet. Secondly, the gravitational torque from the over and 
under-shootings cancel and give almost no net effect. However, note that this is not 
necessarily the case for multi planetary system s or situations in which the satura tion 
of co-rotation torques plays an important role ( Paardekooper &: Papaloizoul |2009| ) . In 



those cases, finite difference codes should only be used with great care and possibly 
with explicit physical viscosity. 

The phenomenon described here does not occur in Godunov codes such as RODEO 
(S.-J. Paardekooper, private communication). This is expected, as those schemes do 
solve the hydro dynamical equations along characteristics and therefore do not have a 
dispersion error. 



E.3.3 Planet torques 

As a further test, we calculate the torques for low mass planets in various disc models 
A planet mass oi q = 1.26 ■ 10~^ and a smoothing length o f b = 0.6 is used, so that the 



simulations can be easily compared to those performed by [Paardekooper fc Papaloizou 



(120091 ). We first test the two dimensional version of Prometheus without physical 
viscosity. The specific torque F = |f x r| felt by the planet is plotted in figure IETtI 
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Figure E.6: Surface density after 10 orbits in simulations with a resolution of 128x384. 
The simulation on the top uses no physical viscosity, whereas the simulation on the 
bottom uses v — 2- 10~^. 
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Figure E.7: Specific torque felt by a planet with mass g = 1.26 ■ 10~^ in an inviscid 
disc with varying surface d e nsity slope a. The results are in good agreement with 
Paardekooper fc Papaloizoij (120091 ) . 



normalised by the square of the planet mass ratio time the disc surface density for four 
different s urface density profiles with slope a. It can be seen by comparing figure IE. 71 to 
figure 7 in lPaardekooper &: Papaloizoul (120091 ) that the measured torques (and therefore 
the migration rates) are in very good agreement. In particular the case a = 3/2 is in 
accord with linear theory (since it does not vary on a libration timescale). The long term 
trend visible towards the end of the simulation is due to different boundary conditions 
used (damping). 



E.3.4 Planet disc interaction in three dimensions 

We now go on and test a three dimensional simulation of pla net disc interactio n. There- 
fore, the results from Prometheus are compared to those of iKley et al.l (120091 ). 

A 20 Earth mass planet is set up on an eccentric and inclined orbit (e = 0.2, i = 5°) 
and embedded in a stratified, three dimensional disc. The disc mass in the computa- 
tional domain (r = 0.4 . . . 2.5) is 0.01 Mq and the aspect ratio of the disc is 0.05. The 
grid resolution used is 128 x 384 x 34. We use no viscosity and a planet potential which 
is smoothed over a fraction (0.43) of a scale height. Note that this is a slightly different 
setup, compared to iKlev et al.l (120091 ). 



Both, the results from Prometheus and those of iKley et al.l (120091 ) are plotted in 
figure IE. 81 They are in good agreement and show only some small variation at large 
times. 
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Figure E.8: Orbital parameters of a planet embedded in a three dimensional disc. The 
solid curves have been obtained wit h Prometheus, the dashed curves have been obtained 
with Nirvana bv lKlev et aD tOQ^ . 
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Figure E.9: Screenshot of the OpenGL visualisation module of Prometheus, showing 
the surface density in a two dimensional simulation with an embedded planet. 

E.4 Visualisations 

Similar to the visualisation module presented in appendix |Dl Prometheus comes with 
a visualisation module, based on the OpenGL API. It can be used to visualise simula- 
tions in real time while they are running or in a standalone playback mode. Standard 
quantities, such as density, velocity, divergence of velocity and vorticity can be plot- 
ted in both, two and three dimensional simulations in either a cylindrical/spherical or 
Cartesian representation. Two screenshots are shown in figures IE. 91 and IE. 101 

The visualisation module is also easily adaptable to allow for direct user interaction 
with the simulation. For example, this can be used to add perturbations or additional 
planets to the simulation and watch the simulation change instantaneously. This pro- 
vides an astonishing new way to easily understand physical phenomena. 

Furthermore, an iPhone version of Prometheus is available on the Apple AppStore 
at http : / / itunes . apple . com/us/app/hydro/ id341643251?mt=8, 
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Figure E.IO: Screenshot of the OpenGL visualisation module of Prometheus, showing 
the density in a three dimensional simulation of a stratified proto-planetary disc with 
an embedded planet. 
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